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Abstract 

Phase  screen  calculations  of  elastic  wave  propagation  in  2-D  random  media  are  com¬ 
pared  with  finite  difference  results  to  assess  the  accuracy  and  efficiency  of  the  former 
method.  The  phase  screen  method  is  a  forward  propagation  algorithm  which  depends 
only  on  the  local  S  and  P  wave  velocities.  It  differs  from  similar  methods  for  scalar 
waves  by  treating  P/S  conversion.  Both  methods  are  used  to  generate  synthetic  seis¬ 
mograms  at  640  evenly  spaced  gridpoints  for  identical  realizations  of  512  x  2750  grids. 
Comparisons  are  made  for  a  suite  of  2-D  random  media  characterized  by  exponen¬ 
tial  and  zeroth  order  von  Karman  (self-similar)  autocorrelation  functions  of  varying 
strength  and  correlation  length.  Constant  and  varying  Poisson  ratio  are  considered. 
Early  arrivals  compare  more  favorably  since  the  phase  screen  method  does  not  include 
4  backscatter.  The  waveforms  are  compared  by  computing  relative  differences  and  cross 

correlations  as  a  function  of  time  offset.  Temporal  energy  centroids  of  bandpass  fil- 
^  tered  synthetic  seismograms  are  also  computed.  Execution  times  are  compared  for 

simulations  on  a  SUN  4/330,  an  ELXSI  6400,  a  CRAY-2,  and  an  nCUBE  parallel 
computer.  The  phase  screen  algorithm  is  roughly  two  orders  of  magnitude  faster  than 
the  finite  difference  algorithm  for  2-D  grids.  A  naive  estimate,  based  on  the  number 
of  operations  in  each  algorithm,  suggests  that  the  phase  screen  method  may  be  three 
orders  of  magnitude  faster  for  3-D  problems. 


1 


1.  INTRODUCTION 


The  purpose  of  this  study  was  to  assess  the  accuracy  and  relative  speed  of 
a  multiple  phase  screen  method  for  propagating  elastic  waves  in  highly  heterogeneous 
media.  The  complexity  of  such  media  makes  analytic  solutions  infeasible.  Thus  we 
have  no  choice  but  to  compare  the  phase  screen  calculations  to  some  other  form  of 
modeled  data.  For  this  comparison,  we  chose  to  use  synthetic  data  calculated  via  the 
finite  difference  method. 

The  finite  difference  technique  is  attractive  because  it  produces  a  full  solu¬ 
tion  to  the  elastic  wave  equation.  Thus  all  direct,  converted,  diffracted  and  guided 
waves  are  accurately  modeled.  Unlike  various  high  frequency  approximations,  there  is 
no  limitation  in  principle  on  the  ratio  of  scatterer  size  to  wavelength.  Also,  synthetic 
seismograms  may  be  generated  at  any  point  on  a  discretized  grid.  The  method  is 
limited  however  by  the  speed  of  the  computer,  the  memory  available  and  the  cost  of 
CPU  time.  With  very  few  exceptions,  the  published  studies  based  on  this  method 
have  been  performed  for  2-D  media  (Frankel,  1989). 

The  phase  screen  method,  developed  by  Fisk  and  McCartor  (1989,1991),  is 
an  algorithm  to  rapidly  forward  propagate  vector  elastic  waves.  Phase  screen  methods 
for  scalar  waves  have  been  used  in  previous  propagation  studies  of  starlight  through 
the  atmosphere  (Ratcliffe,  1956;  Mercier,  1962;  Filice,  1984),  radio  signals  through 
the  ionosphere  (Buckley,  1975;  Bramley,  1977;  Knepp,  1983),  acoustic  waves  in  the 
ocean  (Flatte;  1979),  and  P  waves  in  the  earth  (Haddon  and  Husebye,  1978).  Sim¬ 
ulations  for  scalar  waves  in  three  dimensions  have  been  performed  by  Filice  (1984) 
and  Martin  and  Flatto  (1988).  The  method  used  here  treats  both  S  and  P  waves, 
including  their  conversion,  and  allows  for  simulation  of  seismic  wave  propagation  in 
3-D  heterogeneous  media.  Like  the  finite  difference  method,  this  method  may  be  used 
to  generate  synthetic  seismograms  at  any  gridpoint. 

The  phzise  screen  method  may  be  applied  to  media  described  by  a  mixture 
of  deterministic  and  stochastic  structure.  In  a  previous  study,  Fisk  and  McCartor 
(1991)  compared  the  results  of  the  phase  screen  method  with  an  “exact  solution”  of 
an  elastic  wave  in  a  2-D  laterally-layered  structure.  (By  “exact”  we  mean  that  there 
were  no  fundamental  approximations  made  to  solve  the  equation  of  motion  and  satisfy 
the  boundary  conditions.  The  solution  cannot  be  written  in  closed  form,  however,  but 
may  be  numerically  determined  to  any  desired  precision.)  Even  for  wavelengths  on  the 
order  of  the  length  scale  of  the  layers  and  velocity  variations  of  5%,  the  comparison 
was  excellent. 
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Common  to  ray  theory  methods  (e.g.,  Karal  and  Keller,  1959;  Aki  and 
Richards,  p.  90,  1980;  Ansell,  1981),  the  phase  screen  method  makes  use  of  the  fol¬ 
lowing  facts:  (1)  To  first  order  at  high  frequency  the  P  and  S  wave  displacements, 
expressed  as  the  curl-free  and  divergence-free  portions  of  the  total  displacement  re¬ 
spectively,  decouple;  and  (2)  to  first  order  the  propagation  of  the  waves  depend  only 
on  the  local  P  and  S  wave  velocities,  which  act  to  distort  the  wavefronts.  For  applica¬ 
ble  problems,  the  heterogeneous  medimn  depicted  in  figure  1(a)  may  be  replaced  with 
a  homogeneous  medium  and  a  set  of  phase  screens,  as  shown  in  figure  1(b).  Based 
on  the  first  fact,  the  initial  P  and  S  waves  at  2  =  0  are  independently  propagated 
*  to  z  =  Az  according  to  the  uniform  eleistodynamic  equation,  using  the  average  wave 

velocities.  The  affect  of  variations  in  the  local  velocities  is  trp..ted  by  multiplying  the 
^  S  and  P  waves  by  position  dependent  phase  factors  at  screen  2.  In  general,  the  dis¬ 

torted  wavefronts  no  longer  satisfy  the  curl-free  and  divergence-free  conditions.  Hence, 
if  these  waves  are  decomposed  onto  a  complete  set  of  forward  propagating  P  and  S 
plane  waves  that  do  satisfy  these  conditions,  P/S  conversion  is  obtained.  The  proce¬ 
dure  may  now  be  repeated,  inserting  phase  screens  for  larger  values  of  z,  to  propagate 
the  displacement  further. 

The  spacing  of  the  screens  is  determined  such  that  geometric  optics  applies 
for  propagation  between  them.  Martin  and  Flatt4  (1988)  have  pointed  out  that  diffrac¬ 
tion  effects  are  accumulated  from  propagation  through  many  screens.  By  including 
diffraction  effects,  as  Hudson  (1980)  has  noted,  the  results  are  valid  to  a  much  greater 
range  than  those  of  simple  ray  theory.  The  phase  screen  method  is  limited,  however, 
to  particular  scattering  strengths  and  length  scales  of  the  structure  for  a  given  wave¬ 
length.  Thus  this  method  does  not  have  the  versatility  of  the  finite  difference  method. 
We  believe,  however,  that  it  may  be  used  to  treat  a  significant  range  of  problems  in 
seismology  that  finite  difference  currently  cannot  due  *o  CPU  time  considerations. 

Conversion,  as  Levin  and  Rytov  (1957)  have  noted,  is  a  second  order  effect 
at  high  frequency,  and  is  treated  as  such  in  the  phase  screen  method.  Karal  and  Keller 
(1959),  Richards  (1974)  and  Ansell  (1981)  have  also  pointed  out  that  there  are  second 
order  contributions  to  the  P  and  S  wave  displacements  that  are  not  purely  longitudinal 
'  and  transverse  motion.  These  terms  depend  on  spatial  derivatives  of  the  density  and 

the  Lame  constants.  Thus  treating  P/S  conversion  while  ignoring  these  contributions 
4  is  not  entirely  consistent  to  second  order.  We  have  found,  by  comparison  with  exact 

solutions  for  the  laterally-layered  problem  with  constant  Poisson  ratio,  that  the  results 
are  considerably  more  accurate  by  treating  conversion  in  this  manner  than  ignoring 
it  altogether.  Efforts  to  improve  the  P  and  S  wave  decomposition  are  currently  being 
made. 
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Aside  from  treating  both  P  and  S  wav-}?  in  our  analysis,  our  method  differs 
from  previous  phase  screen  methods  used  in  seisn.  oiogy  by  treating  wide-angle  scatter. 
Although  wide-angle  scatter  has  been  treated  before  in  the  context  of  phase  screen 
methods  for  acoustic  waves  (Thomson  and  Chapman,  1983),  and  for  light  waves  (Feit 
and  Fleck,  1978),  the  typical  starting  point  in  the  formulation  of  the  method  is  to 
assume  the  parabolic  approximation  which  treats  only  small-angle  scatter.  This  ap¬ 
proximation  is  valid  for  either  scalar  or  vector  waves  provided  the  wavelength  is  much 
smaller  than  the  structure  size.  In  a  previous  study,  Fisk  and  McCartor  (1991)  found 
that  the  treatment  of  wide-angle  scatter  allows  for  propagation  in  smaller  structures 
to  be  computed  more  accurately.  Certainly  larger  structures  will  lead  to  better  results 
since  waves  propagating  at  large  oblique  angles  are  apt  to  induce  backscattering.  For 
large  enough  structures,  our  expressions  may  be  expanded  to  recover  the  standard 
parabolic  approximation. 

A  study  by  McLaughlin  and  Anderson  (1987)  hzus  shown  that  methods  that 
do  not  include  wide-angle  scattering  and  conversion,  such  as  simple  ray  theory,  Ry- 
tov  theory  of  a  forward  propagating  scalar  wave,  and  Gaussian  beam  synthesis,  are 
incapable  of  modelling  systematic  delays  of  high  frequencies  (4-5  Hz)  relative  to  lower 
frequencies  (1-2  Hz)  as  has  been  observed  at  NORSAR  in  P-waves  from  Eastern 
Kazakh  explosions.  They  referred  to  this  effect  as  “stochastic  dispersion”  of  the  P- 
wave,  and  quantified  the  time  delay  by  computing  the  temporal  energy  centroids  of 
bandpass  filtered  seismograms.  They  found  that  of  the  methods  they  studied,  only 
finite  difference  calculations  for  media  with  multiple  length  scales  were  able  to  repro¬ 
duce  this  effect.  We  computed  the  temporal  energy  centroids  of  bandpass  filtered  finite 
difference  and  phase  screen  synthetics  to  determine  the  affects  of  including  multiple 
wide-angle  scattering,  but  neglecting  backscatter,  in  the  phase  screen  method. 

The  remainder  of  the  paper  is  organized  as  follows.  Section  2  provides  de¬ 
scriptions  of  the  two  algorithms  used.  Execution  times  are  compared  for  simulations 
on  a  SUN  4/330,  an  ELXSI  6400,  a  CRAY-2,  and  an  nCUBE  parallel  computer.  Due 
to  CPU  time  considerations,  the  finite  difference  simulations  were  only  performed  on 
the  last  two  computers.  Section  3  contains  the  details  of  modeling  the  random  fiuc- 
tuations.  To  directly  compare  the  results  of  the  two  methods,  identical  realizations 
of  the  random  fluctuations  were  used  by  both  routines.  Section  4  contains  compar¬ 
isons  of  the  results.  Synthetics  were  computed  at  640  evenly  distributed  points.  We 
investigate  how  the  accuracy  of  the  phase  screen  method  depends  on  the  autocorre¬ 
lation  function,  the  magnitude  of  the  velocity  fluctuations,  the  correlation  length  of 
the  medium,  and  whether  the  medium  was  described  by  constant  or  varying  Poisson 
ratio.  We  compute  the  cross  correlations  of  the  time  series  and  relative  differences  of 
the  waveforms  for  each  component,  and  average  over  the  receivers  in  the  x-direction. 
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We  also  compare  the  temporal  energy  centroids  of  the  vertical  component  waveforms. 
In  section  5  we  draw  some  conclusions  about  the  applicability  and  usefulness  of  the 
phase  screen  method.  Appendix  A  contains  detailed  formulas  of  the  phase  screen 
method  for  vector  elastic  waves,  and  Appendix  B  addresses  the  accuracy  of  the  finite 
difference  method. 


2.  NUMERICAL  ALGORITHMS 


Finite  Difference  Method 

The  equation  of  motion  for  the  displacement  vector  u  in  a  heterogeneous 
linear  source-free  medium  is 


P(x) 


d^Ui 

'W 


A(x)(V.uK-,+ 


(1) 


where  A(x)  and  /i(x)  are  the  Lame  constants,  and  p(x)  is  the  density.  For  propagation 
in  a  two  dimensions,  the  equations  of  motion  for  the  horizontal  (u)  and  vertical  (v) 
particle  displacements  are 


£ 

dx 


+ 


dz 


,du  dv. 
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dH 

d 

,du  dv. 

d 

.dv  .du 

dx 

{A  +  2m)5j  +  A^ 

(3) 


Throughout  this  work,  we  have  assumed  a  two  dimensional  cartesian  geometry,  and 
that  density  is  constant  throughout  the  medium.  In  addition,  we  have  assumed  the 
medium  is  horizontally  periodic  and  bounded  by  absorbing  boundaries  in  the  vertical 
direction.  (The  latter  boundary  condition  is  irrelevant  to  the  phjise  screen  method, 
since  backscatter  is  ignored.)  The  finite  difference  method  solves  these  coupled  differ¬ 
ential  equations  numerically  by  replacing  the  partial  derivatives  in  space  and  time  by 
finite-difference  approximations.  The  parameters  A,^,p,  and  the  components  of  the 
wavefield  become  functions  of  discrete  2-D  cartesian  coordinates. 

Various  finite  difference  approximations  have  been  used  to  solve  the  wave 
equation  (e.g.,  Boore,  1972;  Alford,  1974;  Kelly  et  al.,  1976;  Virieux,  1986;  Fornberg, 
1987;  Witte,  1989).  For  this  work,  we  have  made  use  of  the  simple  explicit  second- 
order  scheme  proposed  by  Kelly  et  al.  (1976).  Second-order  schemes  are  preferred  in 
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highly  heterogeneous  media  since  they  more  accurately  resolve  rapid  variations  in  the 
wavefield  and  the  medium  (e.g.,  Fornberg,  1987;  Daudt,  1989;  Charrette,  1991).  Also, 
we  favored  a  non-staggered  formulation  because  it  is  easier  to  implement  the  second- 
order  absorbing  boundary  conditions  proposed  by  Clayton  and  Engquist  (1977). 

Phase  Screen  Method 

The  linearity  of  eq.  (1)  allows  the  solution  to  be  written  in  terms  of  P  and 
S  wave  displarcments  as  u  =  Up  +  Us,  which  are  chosen  to  satisfy  the  constraints 
V  X  Up  =  0  and  V  •  U5  =  0.  To  first  order  at  high  frequency  this  decomposition 
results  in  the  decoupling  of  the  equation  of  motion  into  two  simple  wave  equations  of 
the  form 


where  a  time  dependence  of  e”’*"*  is  assumed  for  now.  The  local  propagation  velocities 
for  the  respective  waves  are  given  in  terms  of  the  density  and  Lame  constants  by 
a{x)  =  y'(A(x)  +2Ai(x))/p(x)  and  ^(x)  =  ^A*(x)/p(x). 

The  phase  screen  algorithm  for  vector  waves  may  be  summarized  as  follows. 
(A  more  detailed  discussion  is  piovided  in  Appendix  A.)  Given  a  realization  of  the  2-D 
grid,  it  is  divided  into  equal  intervals  of  length  Az  by  planes  of  constant  z.  Starting 
with  the  initial  displacement,  Up  and  U5  are  uniformly  and  independently  propagated 
to  a  distance  Az.  The  affect  of  the  velocity  perturbations  is  treated  by  multiplying 
the  S  and  P  waves  by  position  dependent  phase  factors,  exp(iA(^))  and  exp(:A^^^). 
The  P-wave  phase  factor  for  the  perturbations  between  z  =  0  and  2  =  Az  is  given  by 
the  geometrical  optics  expression 


A(^)(x,A2)  =  kp  dz  (5) 

Jo  a 

where  kp  =wf  a,  a  is  the  average  P  wave  velocity,  and  6a  is  the  velocity  perturbation. 
There  is  an  analogous  expression  for  the  S-wave  phase  factor.  Similar  expressions  are 
computed  for  the  fiuctuations  between  the  other  screens. 

The  wavefronts,  distorted  by  the  phases,  no  longer  satisfy  the  curl-free  and 
divergence-free  conditions.  Using  standard  Fourier  and  vector  analysis,  these  waves 
are  decomposed  onto  a  complete  set  of  forward  propagating  P  and  S  plane  waves 
that  do  satisfy  these  conditions.  The  procedure  is  repeated,  inserting  phase  screens 
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for  larger  values  of  z,  to  propagate  the  displacement  further.  Synthetic  seismograms 
are  produced  by  computing  the  single-frequency  results  at  a  discrete  set  of  positive 
frequencies.  Using  an  FFT,  the  product  of  single-frequency  results  and  the  Fourier 
spectrum  of  the  time  source  function  is  transformed  to  the  time  domain.  The  real 
part  of  the  transformed  expression  is  then  taken. 

3.  CONSTRUCTING  RANDOM  VELOCITY  PERTURBATIONS 

The  autocorrelation  function  is  commonly  used  to  characterize  random  fields 
and  is  a  measure  for  quantifying  the  similarity  between  neighboring  points  in  a  random 
medium.  It  has  the  property  that  it  is  the  Fourier  transform  of  the  power  spectrum 
<  (Tatarski,  1961).  This  relationship  allows  us  to  build  realizations  from  a  desired 

correlation  function  in  the  wavenumber  domain.  Throughout  this  study,  realizations 
were  constructed  by  convolving  the  square  root  of  the  power  spectrum  with  a  phase 
term  of  the  form  e**,  where  d  is  a  random  number  drawn  from  a  uniform  distribution 
over  the  range  0<$<2n.  Since  the  norm  of  the  phase  term  is  one,  the  shape  of  the 
power  spectrum  and  the  total  power  within  that  spectrum  are  unchanged. 

Three  correlation  functions  have  received  a  great  deal  of  attention  in  the 
scattering  literature;  the  Gaussian,  the  exponential  and  the  von  Karman  functions 
(e.g.  Chernov,  1960;  Tatarski,  1961;  Dainty,  1984;  Frankel  and  Clayton,  1986;  Wu  and 
Aki,  1990).  The  commonly  used  form  of  these  functions  and  their  power  spectra  are 
given  in  Table  1,  and  shown  graphically  in  Figure  2. 

In  the  exponential  function,  the  correlation  length  a  marks  the  lag  where  the 
correlation  function  has  the  value  e~^  (Figure  2).  In  the  wavenumber  domain,  both 
spectra  are  flat  out  to  a  corner  wavenumber  which  is  approximately  equal  to  1/a.  At 
higher  wavenumbers,  the  exponential  falls  off  as  where  N  is  the  number  of 

space  dimensions.  The  fall  off  rate  of  the  spectra  controls  the  amount  of  roughness  in 
the  realization.  Spectra  with  more  energy  at  high  wavenumbers  are  expected  to  show 
more  roughness  (Figure  3)  than  those  which  are  localized  near  zero  wavenumber. 

t  The  von  Karman  function  was  first  introduced  oo  characterize  the  random 

velocity  field  of  a  turbulent  medium  (von  Karman,  1948).  In  the  spatial  domain,  the 
^  von  Karman  function  is  peaked  about  the  origin.  The  peak  is  especially  severe  when 

1/  =  0,  since  then  the  modified  Bessel  functic.^  K^,  goes  to  infinity  as  rf  a  goes  to  zero. 
Although  the  parameter  u  can  take  on  any  value  in  the  range  0  to  1,  it  has  some  special 
properties  at  0,  0.3,  0.5  and  1.  When  1/  =  0  the  spectrum  defines  a  multi-dimensional 
Markov  field  (Goff,  1988),  u  =  0.3  defines  Kolmogorov’s  turbulence  (Wu  and  Aki, 
1990),  while  for  u  —  0.5  the  von  Karman  function  simplifies  to  an  exponential  and 
when  1/  =  1.0  to  an  autoregressive  field. 
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Table  1.  Correlation  Functions  and  Power  Spectra. 


Gaussian  Exponential  von  Karman 


Correlation  Function 


1-D  Power  Spectrum 


2-D  Power  Spectrum 


3-D  Power  Spectrum 


,-r/a 


2a 


l+*?a3 


I* 

T 


(l+fc*o»)s/2 


(a0F)3c-*?“V‘‘ 


1  rrl 
2‘'-»r(i/)  [aj 

“  K^{r/a) 

r(,/  + 1/2) 

rH 

{i+kyy*''^ 

TW  + 1) 

(i 

.+kyy*' 

r(i/  +  3/2) 

8ir’/V 

r(i.) 

(i+kyy*^'^ 

The  peakedness  of  the  von  Karman  correlation  function  leads  to  a  wide 
spectral  representation,  indicating  that  these  media  contain  a  significant  amount  of 
roughness  (Figure  4).  As  in  the  exponential  function,  the  power  spectrum  of  the  von 
Karman  function  is  fiat  up  to  a  corner  wavenumber  roughly  equal  to  1/a.  The  differ¬ 
ence  is  that  at  higher  wavenumbers  the  spectrum  falls  off  as  ^-(^+2*')^  considerably 
slower  than  the  exponential  function.  Thus  for  both  functions,  1/a  defines  a  corner 
wavenumber  and  the  parameter  v  controls  the  rate  of  decay  of  the  power  spectrum 
(Figure  2). 

The  von  Karman  function  has  an  additional  property  that  its  slope  is  dis¬ 
continuous  at  zero  lag.  This  property  qualifies  the  von  Karman  function  as  a  fractal 
(Mandelbrot,  1977).  Fractals  are  unique  and  of  interest  because  they  contain  varia-  * 

tions  on  ail  wavelengths.  Since  many  physical  characteristics  in  the  crust  also  display 
variation  on  a  wide  variety  of  length  scales,  this  autocorrelation  function  may  be  well  ♦ 

suited  to  crustal  applications.  The  self-similar  nature  of  fractals  can  be  easily  seen 
by  examining  the  variance  as  a  function  of  wavenumber.  Figure  5  shows  a  series  of 
1-D  realizations  taken  from  the  three  ACFs  described  above.  All  three  realizations 
have  the  same  correlation  length  (a  =  20  m)  and  were  generated  by  the  same  random 
seed.  At  low  wavenumbers  there  is  little  variation  in  shape  and  variance  between  the 
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traces.  This  is  consistent  with  the  power  spectra  (Figure  2),  which  are  flat  at  low 
wavenumber  for  all  three  functions.  At  high  wavenumbers,  there  is  no  variance  in  the 
Gaussian  trace,  and  the  variance  in  the  exponential  trace  is  smaller  than  it  was  at 
low  wavenumber.  Thus,  for  these  media,  the  variance  over  equal  logzu’ithmic  intervals 
of  wavelength  decreases  as  the  wavelength  decreases  (Frankel,  1989).  This  is  not  so 
for  the  0th  order  von  Karman  function.  The  variance  for  that  function  is  roughly 
constant  over  length  scales  smaller  than  27ra  (Figure  5). 

3.  COMPARISONS 

Details  of  the  Models 

Since  the  most  important  differences  between  the  finite  diflerence  and  phase 
screen  solutions  were  expected  late  in  the  wavetrain,  it  was  necessary  to  simulate 
wave  propagation  over  a  large  enough  region  that  the  scattered  field  could  include 
a  considerable  amount  of  forward  and  backscattered  energy.  This  requirement  led 
to  the  use  of  a  relatively  large  grid  (512  nodes  by  2750  nodes).  A  schematic  of  the 
grid  is  shown  in  Figure  6.  The  spacing  between  the  gridpoints  was  dx  =  dz  =  0.05 
km.  The  lowest  250  gridpoints  in  the  s-direction  were  homogeneous  to  allow  for  the 
pulse  to  be  “initialized”  in  the  finite  difference  code.  The  spatial  form  of  the  source 
was  a  plane  P-wave  located  initially  at  a  depth  of  127.5  km  within  the  homogeneous 
segment.  Arrays  of  64  receivers  were  located  every  250  gridpoints  in  the  s-direction. 
The  phase  screens  were  uniformly  spaced  every  50  gridpoints  in  the  s-direction  for 
the  heterogeneous  portion  of  the  grid. 

The  initial  time  source  function  used  in  our  study  was  a  Ricker  wavelet 
defined  by 

•f’(0  =  (1  -  (27r/o(t  -  to)Y/2)  exp  (-(27r/o(f  -  to))*/^)  .  (6) 

where  /o  is  the  peak  frequency  and  to  defines  the  origin  of  time.  The  value  /o  =  2 
Hz  was  used  for  all  of  the  simulations.  The  finite  difference  simulations  used  10000 
timesteps,  however,  only  1500  of  these  points  were  included  in  the  output.  The  same 
timestep  At  =  0.025s  was  used  for  the  phase  screen  method,  however,  only  256  fre¬ 
quencies  were  sampled;  the  remaining  frequency  data  were  zero  padded  to  minimize 
wrap-around  effects. 

Table  2  lists  the  details  of  the  models,  which  were  chosen  to  explore  how  the 
accuracy  depends  on  the  correlation  function,  the  wavelength  to  correlation  length 
ratio,  and  the  magnitude  of  the  velocity  perturbations.  A  constant  Poisson  ratio 
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was  assumed  for  the  first  four  models.  The  last  four  models  assumed  A  constant; 
thus  the  Poisson  ratio  was  allowed  to  vary  over  the  grid.  The  average  velocities  used 
throughout  this  study  were  a  =  5.5  km/s  and  P  =  afy/S,  Correlation  lengths  of  5 
km  and  2  km  were  used.  The  combination  ka,  relevant  to  the  scattering  regime  (Wu, 
1990),  is  given  in  Table  2,  where  k  is  used  here  to  denote  the  P-wave  wavenumber  at 
peak  frequency.  Although  the  phase  screen  method  is  expected  to  provide  the  most 
favorable  results  for  media  with  only  density  fluctuations,  we  chose  to  consider  only 
variations  in  the  Lame  constants  to  avoid  making  significant  alterations  in  the  finite 
difference  program. 


Table  2.  List  of  Models. 


Model 

ACF 

ka 

1 

Exponential 

2.0% 

2.0% 

11.4 

2 

Self-Similar 

2.0% 

2.0% 

J1.4 

3 

Exponential 

5.0% 

5.0% 

11.4 

4 

Self-Similar 

5.0% 

5.0% 

11.4 

5 

Self-Similar 

1.3% 

2.0% 

11.4 

6 

Self-Similar 

1.3% 

2.0% 

4.6 

7 

Self-Similar 

3.3% 

5.0% 

11.4 

8 

Self-Similar 

3.3% 

5.0% 

4.6 

Comparison  of  Execution  Time 

Although  the  finite  difference  program  had  to  be  modified  to  use  the  nCUBE 
parallel  computer,  the  reduced  execution  time  more  than  compensated  for  the  conver¬ 
sion  of  the  code.  To  take  advantage  of  the  independent  processors  on  the  nCUBE,  the 
large  finite  difference  grid  was  split  into  a  series  of  separate  subgrids,  each  residing  on 
a  separate  nCUBE  processor.  At  each  timestep,  the  independent  subgrids  were  solved 
simultaneously,  the  data  along  the  edges  are  exchanged  and  the  process  is  repeated. 
For  reference,  a  finite  difference  simulation  like  those  presented  here  would  take  more 
than  12  days  to  complete  on  a  fully  dedicated  workstation  (SUN  4/330).  The  same 
problem  could  be  run  on  a  traditional  supercomputer  (CRAY-2)  in  3.2  hrs.  On  a 
128-node  nCUBE  parallel  computer  the  simulation  took  only  1.8  hrs.  to  run,  and  on 
a  1024-node  nCUBE  parallel  computer  the  run  time  was  less  than  1/2  hour. 


10 


Originally  the  phase  screen  code  was  developed  on  an  ELXSI  6400.  This 
version  was  run  on  the  ELXSI  6400,  SUN  4/330,  and  CRAY-2.  To  exploit  the  inde¬ 
pendent  processors  on  the  nCUBE,  minor  modifications  were  made  to  distribute  the 
single-frequency  calculations  among  the  processors.  Table  3  lists  the  run  times  for  the 
phase  screen  simulations.  The  relative  execution  speeds  of  the  phase  screen  and  finite 
difference  simulations  are  roughly  factors  of  125  for  the  SUN  4/330,  52  for  the  CRAY- 
2,  and  63  for  the  128-node  nCUBE  parallel  computer.  (The  phase  screen  code  was 
not  run  on  the  1024-node  nCUBE  machine.)  A  naive  estimate  of  the  relative  speed, 
based  on  the  number  of  operations  in  each  algorithm,  suggests  that  the  phase  screen 
method  may  be  roughly  three  orders  of  magnitude  faster  for  3-D  problems.  Further 
efforts  to  optimize  the  phase  screen  code  on  a  particular  machine  would  likely  lead  to 
even  more  significant  reductions  in  run  times.  Martin  and  Flatte  (1988)  have  shown 
that  vectorizing  a  phase  screen  algorithm  for  scalar  waves  to  run  on  a  CRAY  led  to  a 
speedup  factor  of  nearly  800  as  compared  with  the  run  time  on  a  VAX-11/750. 

Table  3.  Execution  Times  For  The  Phase  Screen  Method. 


Machine 

Number  ol 

Processors 

1 

2 

4 

8 

16 

32 

64 

128 

SUN  4/330 
ELXSI  6400 

8254  s 
1771  s 

saasJ 

CRAY-2 

223  s 

nCUBE 

7964  s 

3954  s 

1997 

s 

1002  s 

510  s 

288  s 

166  s 

103  s 

Comparison  of  Synthetics 

Synthetic  seismograms  for  the  eight  models  are  recorded  for  receivers  located 
at  the  midpoint  in  the  x-direction,  and  for  various  depths  (Figures  7-14).  Figures  7  5-22 
show  comparisons  of  the  synthetics  at  every  sixteenth  receiver  at  a  depth  of  z  =  12.5 
km  [NZ  =  251).  (This  array  of  receivers  was  chosen  since  the  absorbing  boundary 
condition  at  z  =  0  only  partially  absorbs  obliquely  incident  waves  (cf.  Appendix  B). 
Spurious  refiected  signals  are  expected  to  propagate  backward  into  the  grid  in  the  finite 
difference  solution,  however,  they  should  have  smaller  amplitudes  than  the  spurious 
signals  directly  at  the  boundary.)  In  all  of  these  figures  the  solid  and  dotted  curves 
represent  the  finite  difference  and  phase  screen  results,  respectively.  It  is  immediately 
apparent  that  there  is  excellent  agreement  between  the  waveforms  at  early  times, 
which  is  degraded  at  later  times. 
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To  quantify  the  comparison,  relative  differences  in  the  wavefields  were  com¬ 
puted.  The  relative  difference  in  the  peak-to-peak  amplitudes  of  the  first  vertical 
arrivals  were  between  3%  —  6%  for  the  set  of  models.  Figures  23,24  show  the  relative 
differences  in  each  component  of  the  time  series  for  the  eight  models.  We  have  defined 
the  relative  difference  to  be  the  absolute  difference  of  the  phase  screen  and  finite  dif¬ 
ference  traces  at  each  point  in  the  x-direction  and  time  for  z  =  12.5  km,  divided  by 
the  envelope  of  the  finite  difference  amplitudes  for  each  component.  The  results  are 
averaged  over  the  64  receivers  in  the  x-direction  and  smoothed  over  a  running  time 
window  of  T  =  2//o  =  Is.  The  difference  is  particularly  small  for  the  early  arrivals. 
At  later  times  the  waves  become  out  of  phase  with  each  other  and  the  pointwise  com¬ 
parison  worsens.  Because  the  waveforms  oscillate  fairly  rapidly,  even  slight  offsets  in 
phase  can  lead  to  substantial  differences.  The  power  in  the  coda  agrees  well  even 
though  there  are  differences  in  the  relative  phase. 

The  cross  correlation  of  the  wavefields  provides  a  measure  of  the  similarity 
of  the  solutions  as  a  function  of  time  offset.  The  cross  correlation  is  defined  by 
(s(f)s(t  +  r))  for  a  time  series  s(t),  where  r  is  the  time  lag.  The  angled  brackets 
denote  a  statistical  average,  in  this  case  over  the  receivers  in  the  x-direction  and  also 
over  t.  The  cross  correlations  are  normalized  by  the  autocorrelation  at  zero  time 
lag  of  the  finite  difference  results.  The  cross  correlation  for  the  models  are  shown  in 
Figures  25,26.  The  numbers  in  the  figures  adjacent  to  the  model  numbers  correspond 
to  the  ma:wmum  values.  The  wavefields  are  highly  correlated  at  zero  time  lag,  which 
indicates  fEat  the  arrival  times  computed  by  the  two  methods  are  consistent. 

Figures  23-26  provide  a  measure  of  which  types  of  media  the  phase  screen 
method  treats  most  accurately.  Weaker  media  (2%)  are  treated  better  than  the  mod¬ 
erately  strong  media  (5%).  Also,  the  results  compare  more  favorably  for  exponential 
rather  than  self-similar  media.  Backscatter  clearly  plays  a  more  significant  in  the 
latter  cases.  Figure  24  shows  that  the  media  with  a  =  2  km  result  in  only  a  slight 
increase  in  relative  error  as  compared  with  the  media  with  a  =  5  km.  (The  similarity 
in  the  curves  is  a  result  of  using  the  same  random  seed  to  generate  the  media.)  Also, 
the  horizontal  components  are  more  accurate  for  media  with  constant  rather  than 
varying  Poisson  ratio.  This  is  apparently  a  result  of  the  particular  P  and  S  wave 
decomposition  employed  in  the  phase  screen  method. 

Comparison  of  Temporal  Energy  Centroids 

To  measure  the  differences  in  the  treatment  of' conversion  and  wide-angle 
scattering,  including  backscattering,  the  64  synthetics  generated  at  z  =  0  in  each 
model  were  bandpass  filtered  using  a  zero-phase  shift  Gaussian  bandpass  filter  as  in 
the  study  by  McLaughlin  and  Anderson  (1987).  The  filters  were  centered  at  integer 
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frequency  between  1  -  7  Hz  with  full-width  at  half-maximum  of  1  Hz.  The  temporal 
energy  centroids  of  the  unfiltered  and  filtered  synthetics  were  then  computed.  The 
centroid  for  a  time  series  s{t)  is  defined  as 

^  fs(t)Udt 

‘  fs(tydt' 

The  centroids  were  then  averaged  over  the  64  receivers  and  plotted  in  Figures  27-34 
for  the  eight  models.  The  solid  and  dashed  lines  in  the  figures  correspond  to  the 
averaged  centroids  of  the  unfiltered  synthetics.  The  error  bars  represent  ±  one  stan¬ 
dard  deviation  about  the  mean.  For  models  1,2, 5, 6  with  2%  velocity  variations,  the 
comparison  is  quite  good.  There  is  virtually  no  difference  between  the  centroids  of 
the  unfiltered  finite  difference  and  phase  screen  synthetics.  The  time  delay  of  the  4-5 
Hz  frequencies  relative  to  the  1-2  Hz  frequencies  is  exhibited  by  both  methods.  As 
expected,  the  phase  screen  results  do  not  show  quite  as  much  delay,  particularly  for 
models  3,4 ,7,8  where  backscatter  is  more  significant.  Note  that  for  model  3,  where 
the  centroids  do  not  monotonically  increase  with  frequency,  the  phase  screen  cen¬ 
troids  follow  the  same  curve  as  the  finite  difference  centroids,  with  a  slight  offset  due 
to  backscatter.  McLaughlin  and  Anderson  (1987)  attributed  similar  patterns  in  their 
results  to  scattering  resonances.  In  general,  the  amount  of  dispersion  increases  with 
the  magnitude  of  the  velocity  fiuctuations  and  with  the  presence  of  more  small-scale 
structure;  the  self-similar  medium  leads  to  greater  dispersion  than  does  the  exponen¬ 
tial.  There  is  also  an  overall  delay  of  nearly  1  second  of  the  waveforms  in  the  media 
with  5%  fiuctuations  relative  to  the  media  with  2%. 

4.  CONCLUSIONS 

The  random  media  we  have  considered  in  this  study  were  not  intended  to 
model  any  particular  portion  of  the  crust,  however,  they  are  quite  reasonable.  Studies 
of  the  crust  under  the  Montana  LASA  (Aki,  1973;  Capon,  1974;  Bertreussen  et  al., 
1975)  have  used  a  single  layer  isotropic  Gaussian  medium  model.  The  estimated 
thickness  of  the  layer  60  km,  correlation  length  a  =  10  km  and  rms  P  wave  velocity 
perturbation  between  1.9%  —  4%  are  well  within  the  limits  of  our  study.  A  realistic 
model  of  the  crust  and  upper  mantle  under  NORSAR  has  been  proposed  by  Wu  and 
Flatte  (1990).  The  model  consists  of  overlapping  exponential  and  self-similar  layers. 
(Actually,  truncated  power  spectra  corresponding  to  these  correlation  functions  were 
used.)  The  estimated  total  thickness  ~  250  km,  correlation  length  a  >  20  km  and 
rms  P  wave  velocity  perturbation  between  0.5%  —  2.2%  are  also  characteristic  of  our 
study. 
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It  is  evident  that  the  phase  screen  method  is  most  accurate  over  the  first 
portion  of  the  seismograms  where  the  arrivals  have  undergone  only  small-angle  scat¬ 
tering;  at  later  times  backscatter  becomes  more  significant.  Figures  7-22  demonstrate 
the  range  of  accuracy  of  the  phase  screen  method,  depending  on  the  model.  Although 
the  grid  has  been  sufficiently  sampled  such  that  propagation  errors  in  the  finite  differ¬ 
ence  solution  are  negligible,  large-amplitude  spurious  boundary  reflections  are  present 
in  the  finite  difference  synthetics.  Thus,  differences  in  the  solutions  at  intermediate 
to  late  times  can  only  partially  be  attributed  to  errors  in  the  phase  screen  method 
(cf.  Appendix  B). 

We  believe  that  there  are  a  significant  number  of  problems  in  seismology  that  ' 

may  necessitate  the  use  of  the  phzise  screen  method.  It  provides  the  flexibility  to  per¬ 
form  simulations  on  a  variety  of  machines,  from  work  stations  to  supercomputers.  The  j 

substantial  reduction  in  execution  time,  relative  to  the  finite  difference  method,  allows 
for  substantially  longer  propagation  distances,  more  extensive  parametric  studies,  en¬ 
semble  averaging  over  realizations,  and  3-D  simulations  to  be  considered.  Simulations 
in  3-D  media  are  of  particular  interest  if  data  from  three-component  stations,  located 
above  highly  heterogeneous  media,  are  to  be  analyzed. 

Appendix  A.  The  Phase  Screen  Formulation 

Starting  with  the  decoupled  wave  equations  in  expression  (4),  the  solution 
in  three  dimensions,  between  any  consecutive  pair  of  phase  screens,  may  be  expressed 
as 


n{x)  =  /  ^  {8,(kJ  A(kJ  + e'(k^)  B,(kj.)  ,  (A.l) 

where  A  and  Bj  are  Fourier  coefficients  of  the  P  and  S  waves,  respectively.  Adopting 
a  standard  convention  of  summing  over  repeated  indices,  there  is  an  implicit  sum  over 
/  which  labels  the  two  polarizations  of  the  S  wave.  Our  notation  for  the  wave  vectors 
is  as  follows: 


ki 


(/ix)  ^l/)> 


(A.2) 


hp  —  (kx>^a)j  —  (kx,A^)j 


(A.3) 
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ka  =  yjkl  -  k\y  kp  =  \Jk%  -  k]^,  {A.4) 

where  kp  =  |kp|  =  iw/a,  fcg  =  [ksl  =  iw/)3,  and  a  and  ^  are  the  average  propagation 
speeds.  The  unit  vectors  ip  and  ig  are  defined  such  that  the  curl  and  divergence 
constraints  are  satisfied,  i.e.  kp  x  cp  =  0  and  k^  •  ej  =  0,  /  =  1,2.  Unit  vectors 
satisfying  these  conditions  are 


A 

ep 


I 


A  A 

j  X  ks 

\j  X  ks\ 


(A.5) 


(A.6) 


kg  X  eg 


X  cj 


=  {kl  -  ^l)  ^^^ks^{rKky,k\  -  -kykp). 


(A.7) 


Given  an  initial  value  of  the  displacement  at  «  =  0,  the  Fourier  coeffi¬ 
cients  and  Bf^  may  be  determined  by  standard  Fourier  and  vector  projections. 
Once  these  coe  jients  are  known,  eq.  (A.l)  may  be  used  to  propagate  the  displace¬ 
ment  to  z  =  Aj.  Accounting  for  the  local  variations,  the  P  and  S  wave  terms  are 
multiplied  by  phase  factors  of  the  form  in  eq.  (5)  to  obtain 


+  4  k,)  e"'''*  e*'*-+*>*.  (A.8) 


In  a  previous  application  where  the  medium  was  invariant  under  translations  in  the 
a-direction,  the  accumulated  phase  depended  on  the  wavenumber  as  well.  Here  we 
find  that  it  is  adequate  and  more  expedient  to  use  phase  factors  of  the  form  computed 
in  eq.  (5). 

Now  setting  the  initial  value  of  at  z  =  As  equal  to  the  expression  on 
the  RHS  of  eq.  (A.8),  the  new  Fourier  coefficients  and  may  be  determined 
in  terms  of  and  Bf^  and  the  phase  factors.  The  general  form  of  these  relations  is 
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A^^\k,,ky)  =  f  dxdy  [ 

kg  '  kp  •'  ^ 

•  jifcs  •  «p(fcL, ^1)  K) 

=  j-^y-  [  dzdy  f  d^k'  e’(*i“''*)*+‘(*'»-*»)''  (^5  x  e^s{k^tky))  • 

kg  *  kp  *' 

•  |(ibj.  X  Cp(ibi,  Ay)  A(^-^)(Ay  Ay  e‘A‘'*H*.».A^A*) 

+  (fcp  X  «i(AyAy)  flf -'^(AyAy  e’V*  e*A‘®>(*.v.ArA.)|  loj 

The  P/S  coupling  is  evident  in  these  equations.  If  the  phase  factors  are  constant,  it 
is  straightforward  to  show  that  the  waves  decouple.  Thus  conversion  is  caused  by  the 
local  velocities  distorting  the  wavefronts. 

For  the  purposes  of  this  paper,  we  shall  now  consider  only  two  dimensions, 
which  is  recovered  as  a  special  case  by  setting  ky  =  0.  If  Ay  =  0,  the  coefficients 
decouple  from  and  in  the  recursion  relations,  and  hence  will  be  set 
to  zero.  For  computational  purposes,  the  continuous  Fourier  expansions  used  to  ex¬ 
press  the  displacement  are  replaced  with  discrete  versions.  Let  k^  k^  =  27cm/L, 
i4(^)(Ax,Ay)  -»  and  B|^^{A„Ay)  -*  where  m  is  an  integer  and  L  is  the 
length  of  the  screen.  The  displacement  after  N  phase  screens  becomes 

=  X!  A^J^kp,^  +  Y.  (A.ll) 

m  n 

Now  define  the  vector  of  Fourier  coefficients  in  block  form  as 


aL^) 


Bf) 


-  jf")  jl")  «(">  /?(")  B 


(A.12) 
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4 


I 


Mp  and  Ms  are  the  greatest  integers  less  than  wLf2Tra  and  u>L/2jr/5  respectively, 
and  T  denotes  the  transpose  of  the  row  vector.  Only  the  modes  with  harmonic  z- 
dependence  are  included  since  we  are  interested  in  propagation  distances  for  which 
the  evanescent  modes  would  be  exponentially  damped  away.  The  recursion  relations 
for  the  Fourier  coefficients  may  now  be  expressed  as 


(  aL"'  ^ 

Bj")  , 


,  Bf -» 


{A.13) 


The  integers  m,  n  range  from  —Mp  to  +Mp,  and  the  integers  /i,  u  range  from  —Ms 
to  +Mst  and  there  are  implicit  sums  over  riyU  in  this  expression.  In  addition,  the 
A-phases  may  be  Fourier  expanded  as 


g.A('’)(x)  ^  J-  J^(P) 
n 


(A.14) 


and  a  similar  expression  with  P  -*  S.  Combining  this  information  with  eqs.  (A,9,A.10) 
yields 


ttPP  ^fi,rxMa,n  ’I" 

^mn  ““  L-  JL  4.  W 

(A.15) 

TjPS  P  ^mhp,v  —  jj{S) 

»  fc-  jfe  4.  W  ^m-i/ 

(A.16) 

TjSP  ^  ^nha,n  ~  ^nka,n  jj(P) 

(A.17) 

ft 

ttSS  ka,i*kl3,y  +  (5) 

ka  k  +  ik* 

(A.18) 

4 


The  phase  screen  algorithm  for  vector  waves  may  now  be  summarized  as 
follows.  Starting  with  the  initial  displacement  at  z  =  0,  the  initial  vector  of  coeffi¬ 
cients  are  obtained.  Given  a  realization  of  the  2-D  grid,  it  may  be  divided  into  equal 
intervals  of  length  Az  by  planes  of  constant  z.  The  phase  factors  for  each  interval 
are  computed  via  a  discrete  version  of  eq.  (5);  an  inverse  FFT  is  used  to  obtain  the 
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Fourier  coefficients  and  which  are  inserted  into  the  matrix  above.  These 
relations  are  successively  applied  to  the  initial  vector  of  coefficients  until  the  final 
vector  is  obtained.  The  displacement  in  eq.  (A.ll)  is  evaluated  using  an  FFT. 

Appendix  B.  Accuracy  of  the  Finite  Difference  Technique 

The  siccuracy  of  the  finite  difference  formulation  used  here  was  investigated 
by  Charrette  (1991).  In  that  study  it  was  shown  that  the  finite  difference  solution 
converges  to  the  ex2M:t  solution  as  the  discretization  of  the  medium  increases.  In  par¬ 
ticular,  for  the  spatial  and  temporal  sampling  rates  used  in  this  note,  it  was  shown  that 
the  error  in  phase  and  group  velocities  are  less  than  a  few  percent.  The  accuracy  tests 
described  in  Charrette  (1991)  assume  the  medium  is  largely  homogeneous.  When  the 
medium  is  highly  heterogeneous,  errors  in  calculating  the  velocity  gradients  will  also 
effect  the  solution.  The  behavior  of  the  finite  difference  technique  is  these  situations  is 
not  well  understood  and  is  likely  to  depend  greatly  on  the  nature  of  the  medium.  Due 
to  the  complexity  of  the  problem,  the  only  tractable  way  to  estimate  these  errors  is 
to  exploit  the  fact  that  in  the  co.'tinuum  limit  the  finite  difference  solution  converges 
to  the  true  solution. 

To  illustrate  this,  we  ran  one  simulation  at  twice  the  spatial  and  temporal 
sampling  rate  used  throughout  this  study.  This  resulted  in  a  sampling  rate  of  55  points 
per  wavelength  for  the  shortest  wavelength  P-wave  (110  points  per  wavelength  for  the 
center  frequency)  and  a  sampling  rate  of  30  points  per  wavelength  for  the  shortest 
wavelength  S-wave  (60  points  per  wavelength  for  the  center  frequency).  Figures  31 
and  32  show  that  the  difference  between  the  two  solutions  is  negligible. 

Additional  errors  are  induced  into  the  finite  difference  solution  from  the 
boundary  conditions  at  the  ends  of  the  grid.  Randall  (1988)  has  noted  that  parax¬ 
ial  boundary  conditions,  such  as  the  one  proposed  by  Clayton  and  Engquist  (1977), 
perfectly  absorb  elastic  waves  normally  incident,  but  only  partially  absorb  obliquely 
incident  waves.  He  further  stated  that  portions  of  the  wave  incident  at  angles  between 
about  30°  and  45°  may  lead  to  unacceptable  boundary  reflections.  The  horizontal  dis¬ 
placement  is  apt  to  contain  more  large-amplitude  spurious  reflections  (e.g..  Fig.  10, 
Frankel  and  Clayton,  1986;  Fig.  6,  Randall,  1988).  Also,  stronger  scattering  media  are 
more  likely  to  produce  spurious  signals  since  there  will  typically  be  more  wide-angle 
scatter  incident  on  the  boundary. 

ACKNOWLEDGEMENTS 

One  of  the  authors  (EEC)  would  like  to  thank  his  employer,  nCUBE  Corporation, 
for  permission  to  participate  in  this  publication.  Many  of  the  computer  simulations 


18 


presented  in  this  study  were  performed  on  the  192-node  nCUBE2  computer  at  the 
Geophysical  Center  for  Parallel  Processing,  part  of  the  Earth  Resources  Laboratory 
at  MIT.  We  would  like  to  thank  the  members  of  that  lab  for  allowing  us  to  use  a 
considerable  amount  of  disk  space.  This  work  was  supported  by  the  Defense  Advanced 
Research  Projects  Agency  through  contract  F19628-89-C-0040  administered  by  the  Air 
Force  Geophysical  Laboratory. 


REFERENCES 

Aki,  K.,  Scattering  of  P  waves  under  the  Montana  Lasa,  J.  Gtophys.  Rea.,  78,  1334- 
1246,  1973. 

Aki,  K.,  and  P.  G.  Richards,  Quantitative  Seismology:  Theory  and  Methods,  Vols.  1 
and  2,  W.  H.  Freeman,  San  Francisco,  Calif.,  1980. 

Alford,  R.  M.,  K.  R.  Kelly,  and  D.  M.  Boore,  Accuracy  of  finite  difference  modeling 
of  the  acoustic  wave  equation.  Geophysics,  39,  834-842,  1974. 

Ansell,  J.  H.,  The  decoupling  of  P  and  S  waves  in  inhomogeneous  media.  Physics  of 
the  Earth  and  Planetary  Interiors,  26,  246-249,  1981. 

Bertreussen,  K.  A.,  A.  Christoffersson,  E.  S.  Husebye,  and  A.  Dahle,  Wave  scattering 
theory  in  analysis  of  P  wave  anomalies  at  NORSAR  and  LASA,  Geophya.  J.  R.  Aatr.  Soc., 
42,  402-417,  1975. 

Boore,  D.  M.,  Finite  difference  methods  for  seismic  wave  propagation  in  heterogeneous 
materials,  pp.  1-37,  vol,  11  of  Methods  in  Computational  Physics,  Academic  Press, 
New  York,  1972. 

Bramley,  E.  N.,  The  accuracy  of  computing  ionospheric  radio-wave  scintillation  by 
the  thin-phase  screen  approximation,  J.  Atmos.  Terr.  Phys.,  39,  367-373,  1977. 

Buckley,  R.,  Diffraction  by  a  random  phase-changing  screen:  A  numerical  experiment, 

J.  Atmos.  Terr.  Phys.,  37,  1431-1446,  1975. 

Capon,  J.,  Characterization  of  crust  and  upper  mantle  structure  under  LASA  as  a 
random  medium.  Bull.  Seismol.  Soc.  Am.,  64,  235-266,  1974. 

Charrette,  E.  E.,  Elastic  Wave  Scattering  in  Random  Media,  Ph.D.  Thesis,  Mas¬ 
sachusetts  Institute  of  Technology,  Cambridge,  MA,  1991. 

Chernov,  L.  A.,  Wave  Propagation  in  a  Random  Medium,  McGraw-Hill,  New  York, 
1960. 


19 


Clayton,  R.  and  B.  Engquist,  Absorbing  boundary  conditions  for  acoustic  and  elastic 
wave  equations,  Bull.  Seismol.  Soc.  Am.,  67,  1529-1540,  1977. 

Dainty,  A.  M.,  High-frequency  acoustic  backscattering  and  seismic  attenuation,  J.  Gto- 
pkys.  Res.,  89,  3172-3176, 1984. 

Daudt,  C.  R.,  L.  W.  Braile,  R.  L.  Nowack,  and  C.  S.  Chiang,  A  comparison  of  finite- 
difiFei'ence  methods.  Bull.  Seismol.  Soe.  Am.,  79,  1210-1230,  1989. 

Feit,  M.  D.  and  J.  A.  Fleck,  Jr.,  Light  propagation  in  graded-index  optical  fibers. 
Applied  Optics,  vol.  17,  no.  24,  3990-3998,  1978. 

Filice,  J.  P.,  Studies  of  the  Microscale  Density  Fluctuations  in  the  Solar  Wind  using 
the  Power  Law  Phase  Screen  Model,  U.  Calif.,  San  Diego,  Ph.D.  Thesis,  1984. 

Fisk,  M.  D.,  and  G.  D.  McCartor,  “The  Phase-Screen  Method  for  Elastic  Waves  and 
Seismic  Discrimination”,  GL-TR-89-0330,  ADA220771,  Mission  Research  Corp.,  1989. 

Fisk,  M.  D.,  and  G.  D.  McCartor,  The  phase-screen  method  for  vector  elastic  waves, 
J.  Geophys.  Res.,  96,  5985-6010, 1991. 

Flatt6,  S.  M.  (Editor),  Sound  Transmission  Through  A  Fluctuating  Ocean,  Cambridge 
University  Press,  1979. 

Fornberg,  B.,  The  pseudospectral  method:  Comparisons  with  finite  differences  for  the 
elastic  wave  equation.  Geophysics,  52,  483-501,  1987. 

Frankel,  A.,  A  review  of  numerical  experiments  on  seismic  wave  scattering.  Pure 
Appl.  Geophys.,  131,  639-685,  1989. 

Frankel,  A.,  and  R.  W.  Clayton,  Finite  difference  simulations  of  seismic  scattering: 
Implications  for  the  propagation  of  short-period  seismic  waves  in  the  crust  and  models 
of  crustal  heterogeneity,  J.  Geophys  Res.,  91,  6465-6489,  1986. 

Goff,  J.  A.  and  T.  H.  Jordan,  Stochastic  modeling  of  seafloor  morphology:  inversion 
of  sea  beam  data  for  second  order  statistics,  J.  Geophys.  Res.,  9S,  13589-13608, 1988. 

Haddon,  R.  A.  W.  and  E.  S.  Husebye,  Joint  interpretation  of  P-wave  time  and  am¬ 
plitude  anomalies  in  terms  of  lithospheric  heterogeneities,  Geophys.  J.  R.  Astr.  Soc., 
55,  19-43,  1978. 

Hudson,  J.  A.,  A  parabolic  approximation  for  elastic  waves.  Wave  Motion,  2,  207-214, 
1980. 


20 


Karal,  F.  C.,  and  J.  B.  Keller,  Elastic  wave  propagation  in  homogeneous  and  inhomo¬ 
geneous  media,  J.  Acoust.  Soc.  Am..,  31,  694-705,  1959. 

Kelly,  K.  R.,  R.  W.  Ward,  S.  Treitel,  and  R.  M.  Alford,  Synthetic  seismograms:  A 
finite  difference  approach.  Geophysics,  41,  2-27,  1976. 

Knepp,  D.  L.,  Multiple  phase  screen  calculation  of  the  temporal  behavior  of  stochastic 
waves,  Proc.  IEEE,  71,  722-737,  1983. 

Levin,  M.  K.,  and  S.  M.  Rytov,  Transitions  to  geometrical  approximations  in  elzisticity 
theory,  Sov.  Phys.  Acoust.,  2,  179-184,  1957. 

Mandelbrot,  B.  B.,  Fractals,  W.  H.  Freeman,  San  Francisco,  1977. 

McLaughlin,  K.  L.,  and  L.  M.  Anderson,  Stochastic  dispersion  of  short-period  P-waves 
due  to  scattering  and  multipathing,  Geophys.  J.  R.  Astr.  Soc.,  89,  933-963,  1987. 

Martin,  J.  M,  and  S.  M.  Flatte,  Intensity  images  and  statistics  from  niunerical  sim¬ 
ulation  of  wave  propagation  in  3-D  random  media.  Applied  Optics,  27,  2111-2125, 
1988. 

Mercier,  R.  P.,  Diffraction  by  a  screen  causing  large  random  phase  fluctuations, 
Proc.  Cambridge  Phil.  Soc.,  58,  382-400,  1962. 

Randall,  C.  J.,  Absorbing  boundary  conditions  for  the  elastic  wave  equation.  Geo¬ 
physics,  53,  611-624,  1988. 

Ratcliffe,  J.  A.,  Reports  on  Progre."'!  in  Physics,  XIX,  190-263,  1956. 

Richards,  P.  G.,  Weakly  coupled  potentials  for  high-frequency  elastic  waves  in  contin¬ 
uously  stratified  media.  Bull.  Seismol.  Soc.  Am.,  64,  1575-1588,  1974. 

Tatarski,  V.  I.,  Wave  Piopagation  in  a  Turbulent  Medium,  McGraw-Hill,  New  York, 
1961. 

Thomson,  D.  J.  and  N.  R.  Chapman,  A  wide-angle  split-step  algorithm  for  the  parabolic 
equation,  J.  Acoust.  Soc.  Am.,  74,  1848-1854,  1983. 

Virieux,  J.,  P-SV  wave  propagation  in  heterogeneous  media:  velocity-stress  finite- 
difference  method.  Geophysics,  51,  889-901,  1986. 

von  Karman,  T.,  Progress  in  the  statistical  theory  of  turbulence,  J.  Mar.  Res.,  7, 
252-264,  1948. 


21 


Witte,  D.,  The  Pseudospectral  Method  for  Simulating  Wave  Propagation,  Ph.D.  The¬ 
sis,  Columbia  University,  Palisades,  NY,  1989. 

Wu,  R.  S.  and  K.  Aki,  The  perturbation  method  in  elastic  wave  scattering.  Pure 
Appl.  Geophys.,  ISl,  605-637, 1990. 

Wu,  R.  S.  and  S.  M.  Flatt6,  Transmission  fluctuations  across  an  array  and  hetero¬ 
geneities  in  the  crust  and  upper  mantle,  Pure  Appl.  Geophys.,  1S2,  175-196,  1990. 

Wu,  R.  S.,  Seismic  wave  scattering,  in  Encyclopedia  of  Geophysics,  D.  E.  James, 

ed.  Van  Nostrand  Reinhold,  1988.  * 


# 


22 


z  =  0 


2  =  Az 


Figure  1.  The  multiple  phase-screen  method  replaces  (a)  each  segment 
of  the  heterogeneous  medium  with  (b)  a  uniform  segment  and 
a  phase  screen.  The  accumulated  position-dependent  phase  is 
projected  onto  screen  2. 
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The  model  autocovariance  functions  (top)  and  their  1>D 
power  spectra  (bottom).  The  spectra  and  normalized  so 
that  they  have  the  same  power. 
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Figure  3.  A  2-D  realization  of  a  random  medium  with  an  exponen¬ 
tial  autocorrelation  function.  The  correlation  length  in  this 
realization  is  20m,  and  there  is  5%  RMS  deviation  in  the 
velocity. 
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Figure  4. 


Same  as  Figure  2, 1;  with  a  0th  order  von  Karman  auto> 
correlation  function  ( Jharrette,  1991). 
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Figure  5.  Random  realizations  from  the  1-D  Gaussian,  exponential,  and  0th  order  von  Karman 
autocorrelation  functions,  a)  unfiltered,  b)  bandpass  filtered  allowing  wavelengths 
2.5a-5a,  c)  bandpass  filtered  allowing  wavelengths,  a/4>a/2.  All  realizations  were 
constructed  with  the  same  random  seed,  and  are  plotted  at  constant  scale  (After 
(Frankel,  1989)). 
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Schematic  of  the  2-D  grid.  Arrays  of  64  receivers  were 
located  every  250  gridpoints  in  the  2-direction. 
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Figure  7.  Model  1  gynthetics  of  the  vertical  (V)  and  horisontal  (U) 
components  of  the  displacement.  Receiver  locations  were 
at  the  x-direction  midpoint  for  three  depths. 
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Figure  8.  Model  2  synthetics  of  the  vertical  (V)  and  horizontal  (U) 
components  of  the  displacement.  Receiver  locations  were 
at  the  z-direction  midpoint  for  three  depths. 
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Figure  10.  Model  4  synthetics  of  the  vertical  (V)  and  horisontal  (U) 
components  of  the  displacement.  Receiver  locations  were 
at  the  ^-direction  midpoint  for  three  depths. 


32 


MODEL  #5 


COORD 

V 

257,1 

U 

V 

257,251 

U 

V 

257,501 

U 


15  20  25  30  35 

t  (sec) 


Figure  11,  Model  S  synthetics  of  the  vertical  (V)  and  horizontal  (U) 
components  of  the  displacement.  Receiver  locations  were 
at  the  x-direction  midpoint  for  three  depths. 
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Figure  12.  Model  6  synthetics  of  the  vertical  (V)  and  horizontal  (U) 
components  of  the  displacement.  Receiver  locations  were 
at  the  x-direction  midpoint  for  three  depths. 
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Figure  13.  Model  7  synthetics  of  the  vertical  (V)  and  horizontal  (U) 
components  of  the  displacement.  Receiver  locations  were 
at  the  x-direction  midpoint  for  three  depths. 
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Figure  14.  Model  8  synthetics  of  the  vertical  (V)  and  horizontal  (U) 
components  of  the  displacement.  Receiver  locations  were 
at  the  z-direction  midpoint  for  three  depths. 
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Fig.'  e  15.  Model  1  sythetics  for  4  evenly  spaced  receivers  at  a  depth  of  12.5  km. 
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Figure  16.  Model  2  sythetics  for  4  evenly  spaced  receivers  at  a  depth  of  12.5  km. 
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Figure  18.  Model  4  sythetics  for  4  evenly  spaced  receivers  at  a  depth  of  12.5  km. 
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Figure  19.  Model  5  sythetics  for  4  evenly  spaced  receivers  at  a  depth  of  12.5  km. 
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Figure  20.  Model  6  sythetics  for  4  evenly  spaced  receivers  at  a  depth  of  12.5  km. 
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Figure  21.  Model  7  sythetics  for  4  evenly  spaced  receivers  at  a  depth  of  12.5  km. 
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Figure  22.  Model  8  sythetics  for  4  evenly  spaced  receivers  at  a  depth  of  12.5  km. 
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Figure  24.  Relative  differences  in  the  phase  screen  and  finite  difference 
synthetics  for  models  5-8. 
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Figure  25. 
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Cross  correlations  as  a  function  of  time  lag  between  the 
phase  screen  and  finite  difference  synthetics  for  models  1>4. 
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Figure  26. 
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Cross  correlations  as  a  function  of  time  lag  between  the 
phase  screen  and  finite  difference  synthetics  for  models  5-8. 
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TEMPORAL  ENERGY  CENTROIL  J:  MODEL  #1 


Figure  27.  Average  temporal  energy  centroids  of  bandpass  filtered  synthetics  at  z  =  0  for  model 
1.  The  solid  and  dashed  lines  correspond  to  the  averaged  centroids  of  the  unfiltered 
finite  difference  and  phase  screen  synthetics,  respectively. 
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Figure  28.  Average  temporal  energy  centroids  of  bandpass  filtered  sy  nthetics  at  2  =  0  for 
model  2. 
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Figure  29.  Average  temporal  energy  centroids  of  bandpass  filtered  synthetics  at 
model  3. 
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Figure  30.  Average  temporal  energy  centroids  of  bandpass  filtered  synthetics 
model  4. 
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Figure  31.  Average  temporal  energy  centroids  of  bandpass  filtered  synthetics  at 
model  5. 
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Figure  32.  Average  temporal  energy  centroids  of  bandpass  filtered  synthetics  at  z  —  0  for 
model  6. 
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Figure  33.  Average  temporal  energy  centroids  of  bandpass  filtered  synthetics  at 
model  7. 
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Figure  34.  Average  temporal  energy  centroids  of  bandpass  filtered  synthetics  at 
model  8. 
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Figure  36.  Comparison  of  synthetic  seismograms  using  two  different  grid  spacings.  Shown  here 
is  the  horizontal  component  of  the  displacement  vector  (dominated  by  shear  waves) 
for  model  1.  The  solid  lines  are  from  the  finite  difference  simulations  with  50  m  grid 
spacings  and  the  dashed  lines  are  from  the  finite  difference  simulations  with  25  m 
grid  spacings. 
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Abstract 

Phase  screen  simulations  of  vector  wave  propagation  in  elastic  random  media  are 
applied  to  two  independent  studies  of  seismic  wave  scattering.  In  the  first  study, 
we  analyze  transmission  fluctuations  of  incident  P  waves  propagated  through  a  two- 
dimensional  version  of  Flatt^  and  Wu’s  two-layer  model  of  the  crust  and  upper  mantle 
under  NORSAR.  Transverse  coherence  functions  (TCF’s)  of  the  log  amplitude,  phase 
and  cross  correlation,  computed  from  vector  wave  simulations  using  100  realizations, 
are  compared  to  analytic  and  simulated  TCF’s  based  on  the  parabolic  approximation 
to  the  scalar  wave  equation;  the  analytic  results  also  assume  the  Rytov  approximation. 
Our  results  show  that  the  TCF’s  all  agree  reasonably  well  for  velocity  perturbations 
<  0.5%,  but  the  simulated  TCF’s  (vector  and  scalar)  are  significantly  different  for 
perturbations  >  2%.  The  simulated  vector  wave  also  decorrelates  over  somewhat 
shorter  spatial  offsets  than  the  simulated  scalar  wave  for  perturbations  above  2%. 
In  the  second  study,  we  examine  the  effect  of  random  structure  on  the  variances  of 
the  direct  peak-to-peak  P  wave  amplitude,  and  the  rms  amplitude  of  the  transverse 
component  in  the  velocity  window  between  the  P  and  S  wave  speeds.  A  large  synthetic 
database  of  seismograms  were  generated  at  25  km  intervals  out  to  200  km  for  each 
of  the  eight  random  medixim  models  considered.  Our  study  shows  that  the  relative 
variance  of  the  scattered  phase  is  far  less  dependent  on  the  structure  than  the  direct 
ph£ise.  The  variance  of  the  forward  scattered  P  wave  depends  greatly  on  the  strength 
of  the  large-scale  heterogeneities. 


1.  INTRODUCTION 

The  complexity  of  the  heterogeneities  in  the  lithosphere  and  asthenosphere 
has  led  to  a  statistical  characterization  of  the  small-scale  structure  as  random  me¬ 
dia.  Wave  propagation  in  random  media  is  a  relatively  mature  subject  (e.g.,  Chernov, 
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1960;  Tatarskii,  1961;  Ishimaru,  1977;  Flatt6  et  al.,  1979).  Predominantly  the  theory 
has  been  applied  to  scalar  wave  propagation.  Typically,  weak  scattering  approxi¬ 
mations  are  assumed  to  make  the  calculations  tractable.  For  propagation  in  highly 
heterogeneous  media,  however,  analytic  solutions  are  often  infeasible.  Numerical  sim¬ 
ulations,  such  as  finite  difference  calculations,  are  attractive  in  this  case  because  they 
produce  complete  solutions  to  the  elastic  wave  equation,  and  synthetic  data  may  be 
generated  at  any  gridpoint  in  the  model  (Frankel,  1989).  Unfortunately,  CPU  time 
considerations  often  prohibit  the  use  of  finite  difference  techniques  for  many  interesting 
problems  in  seismology.  For  the  case  of  moderately  strong  scattering,  a  practical  way 
to  obtain  the  solution  is  to  make  a  compromise  between  the  analytical  and  numerical 
methods,  e.g.,  by  neglecting  backscatter,  numerical  simulations  may  be  performed  in 
a  fraction  of  the  time  required  by  full  finite  difference  calculations. 

In  this  paper,  we  consider  two  independent  problems  of  seismic  wave  scat¬ 
tering  in  random  media,  whose  analytic  and  finite  difference  solutions  are  currently 
infeasible.  This  is  not  to  say  that  methods  of  these  types  could  not  be  used  if  certain 
assumptions  were  made  about  the  strength  of  the  scattering,  or  the  extent  of  the 
medium.  However,  the  interesting  features  of  these  problems  cannot  be  adequately 
probed  with  these  solutions.  In  the  first  study,  we  analyze  transmission  fluctuations 
of  incident  plane  P  waves,  propagated  through  a  two-dimensional  version  of  the  model 
proposed  by  Flatt4  and  Wu  (1988)  for  the  crust  and  upper  mantle  beneath  NORSAR. 
Analytic  solutions  exist  for  this  problem,  however,  the  approximations  involved  lead  to 
results  that  are  valid  only  for  weak  fluctuations.  We  compare  simulated  and  analytic 
solutions  to  assess  the  validity  of  the  approximations  assumed,  and  the  usefulness  of 
phase  screen  simulations  for  analyzing  array  data.  In  the  second  study,  we  examine 
the  effect  of  random  structure  on  the  variances  of  direct  and  scattered  phases.  A 
large  statistical  ensemble  of  realizations  (320  total  samples  of  synthetic  seismograms 
at  25  km  intervals  out  to  200  km  for  each  of  eight  models)  was  generated  to  compute 
statistically  reliable  variances.  It  is  the  first  study  of  this  type  of  which  we  are  aware. 
The  details  of  these  studies  are  contained  in  sections  2  and  3.  First,  we  briefly  describe 
the  computational  technique  used  for  this  work. 

1.1.  The  Phase  Screen  Method 

The  phase  screen  method  is  a  rapid  forward  propagation  algorithm,  which 
exploits  the  benefits  of  analytic  approximations  and  numerical  simulation.  Phase 
screen  methods  for  scalar  waves  have  been  used  in  previous  propagation  studies  of 
starlight  through  the  atmosphere  (Ratcliffe,  1956;  Mercier,  1962;  Filice,  1984),  radio 
signals  through  the  ionosphere  (Buckley,  1975;  Bramley,  1977;  Knepp,  1983),  acoustic 
waves  in  the  ocean  (Flatte  et  al.,  1979),  and  P  waves  in  the  crust  (Haddon  and  Huse- 
bye,  1978).  Simulations  for  scalar  waves  in  three  dimensions  have  been  performed  by 
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Filice  (1984)  and  Martin  and  Flatty  (1988).  Recently,  Fisk  and  McCartor  (1989,1991) 
have  developed  a  phase  screen  method  which  treats  vector  elastic  waves,  including 
P/S  conversion.  The  algorithm  allows  for  simulation  of  seismic  wave  propagation  in 
3-D  heterogeneous  media.  Like  the  finite  difference  method,  this  method  may  be  used 
to  generate  synthetic  data  at  any  gridpoint. 

Common  to  ray  theory  methods  (e.g.,  Aki  and  Richards,  1980),  the  phase 
screen  method  makes  use  of  the  following  facts:  (1)  To  first  order  at  high  frequency  the 
P  and  S  wave  displacements,  expressed  as  the  curl-free  and  divergence-free  portions 
of  the  total  displacement  respectively,  decouple;  and  (2)  to  first  order  the  propagation 
of  the  waves  depend  only  on  the  local  P  and  S  wave  velocities,  which  act  to  retard 
or  advance  the  phases  of  the  wavefronts.  For  problems  which  meet  these  criteria,  the 
heterogeneous  medimn  depicted  in  Figure  1(a)  may  be  replaced  with  a  homogeneous 
medium  and  a  set  of  phase  screens,  as  shown  in  Figure  1(b).  Based  on  the  first 
fact,  and  working  with  single-frequency  plane  waves  for  now,  the  initial  P  and  S 
waves  at  2  =  0  are  independently  propagated  to  z  =  Az  according  to  the  uniform 
elastodynamic  equation,  using  the  average  wave  velocities.  The  affect  of  variations  in 
the  local  velocities  is  treated  by  multiplying  the  S  and  P  waves  by  position  dependent 
phase  factors  at  screen  2.  The  spacing  of  the  screens  is  determined  such  that  the  phase 
factors  may  be  computed  from  geometrical  optics.  Below,  we  describe  how  random 
realizations  of  the  phase  factors  are  computed. 

The  wavefronts  of  the  P  and  S  waves  become  distorted  by  the  phase  factors, 
and  thus  no  longer  satisfy  the  curl-free  and  divergence-free  conditions.  Hence,  if  these 
waves  are  decomposed  onto  a  complete  set  of  forward  propagating  P  and  S  plane 
waves  that  do  satisfy  these  conditions,  P/S  conversion  occurs.  The  procedure  may 
now  be  repeated,  inserting  phase  screens  for  larger  values  of  z,  to  propagate  the  new 
expressions  for  the  P  and  S  wave  displacements  further.  If  synthetic  data  in  the  time 
domain  are  desired,  the  single-frequency  results  may  be  computed  at  a  discrete  set 
of  frequencies,  and  then  convolved  with  the  frequency  spectrum  of  the  time  source 
function.  (Detailed  formulas  of  the  method  outlined  here  have  been  given  by  Fisk  and 
McCartor  (1989,1991)  and  Fisk  et  al.  (1991).) 

The  method  has  been  compared,  with  considerable  success,  to  other  calcu¬ 
lations.  Fisk  and  McCartor  (1991)  compared  the  results  of  the  phase  screen  method 
to  an  “exact  solution”  of  an  elastic  wave  in  a  2-D  laterally-layered  structure.  (By 
“exact”  we  mean  that  there  were  no  fundamental  approximations  made  to  solve  the 
equation  of  motion  and  satisfy  the  boundary  conditions.  The  solution  could  not  be 
written  in  closed  form,  however,  but  may  be  numerically  determined  to  any  desired 
precision.)  Even  for  wavelengths  on  the  order  of  the  length  scale  of  the  layers  and 
velocity  variations  of  5%,  the  comparison  was  excellent.  In  a  separate  study,  Fisk  et 
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al.  (1991)  compared  phase  screen  and  finite  difference  calculations  for  elastic  waves  in 
2-D  random  media.  Again,  the  comparison  was  quite  favorable.  In  addition,  it  was 
shown  that  the  phase  screen  algorithm  is  roughly  two  orders  of  magnitude  faster  than 
a  second-order  finite  difference  scheme  for  two-dimensional  problems. 

2.  TRANSMISSION  FLUCTUATIONS 

Transmission  fluctuations  of  direct  P  arrivals  have  been  analyzed  at  LASA 
and  NORSAR  by  Aki  (1973),  Capon  (1974),  Bertreussen  et  al.  (1975)  and  others  to 
deduce  the  statistical  properties  of  the  crust  and  upper  mantle  beneath  the  arrays. 
Assuming  a  single-layer  Gaussian  medium,  comparisons  of  theoretical  TCP’s  of  the 
log  amplitude  and  phase  with  data,  provided  estimates  of  the  depth,  correlation  length 
and  rms  velocity  perturbations  of  the  crust.  The  TCP’s  measure  how  the  waveform 
decorrelates  as  a  function  of  the  spatial  offset  across  an  array,  due  to  propagation 
through  the  random  structure. 

The  analytic  techniques  to  compute  TCP’s  were  originally  developed  by 
Chernov  (1960),  Tatarskii  (1971),  Munk  and  Zachariasen  (1976),  Ishimaru  (1978), 
and  Plattd  et  al.  (1979).  Platte  and  Wu  (1988)  extended  the  theory  to  compute  an¬ 
gular  coherence  functions  (ACP’s),  which  describe  the  coherence  of  the  wavefield  as  a 
function  of  the  angle  of  incidence.  Using  the  additional  information  in  the  ACP’s,  they 
found  that  the  single-layer  Gaussian  medium  model  does  not  accurately  fit  the  data 
at  NORSAR.  They  proposed  a  two-layer  power-law  medium  model  which  sufficiently 
predicts  the  observations. 

While  the  analytic  calculations  provide  elegant  closed  form  expressions  for 
the  coherence,  they  often  assume  the  scalar  wave,  Rytov  and  parabolic  approxima¬ 
tions.  Wu  and  Platte  (1990)  have  noted  that  their  analytic  results  are  valid  for  only 
sufficiently  weak  fluctuations.  In  fact,  by  assuming  the  Rytov  approximation,  the  nor¬ 
malized  TCP’s  are  independent  of  the  magnitude  of  the  rms  velocity  perturbations. 
Intuitively,  the  coherence  of  the  wave  cannot  be  entirely  independent  of  the  scattering 
strength  of  the  medium.  Wu  and  Platt4  (1990)  suggested  that  numerical  simulation 
is  needed  to  determine  when  the  weak  fluctuation  (Rytov)  approximation  no  longer 
applies. 


The  Rytov  approximation  assumes  that  gradients  in  the  wave  fluctuations 
are  small.  Thus,  as  the  frequency  of  the  wave  increases,  the  validity  of  the  approxima¬ 
tion  decreases.  Using  a  phase  screen  method  based  on  the  parabolic  approximation 
to  the  scalar  wave  equation,  Platte  et  al.  (1988)  looked  at  the  behavior  of  the  TCP’s 
with  increasing  frequency.  They  found  that  at  2  Hz,  their  simulated  and  analytic  re¬ 
sults  were  in  agreement.  Por  frequencies  greater  than  4  Hz,  the  log  amplitude  TCP’s 
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differed  dramatically.  Alternatively,  as  the  strength  of  the  perturbations  increase,  for 
fixed  frequency,  the  analytic  results  are  also  expected  to  break  down.  The  scalar  wave 
and  parabolic  approximations  are  also  apt  to  break  down  in  this  case  l  nee  moder¬ 
ately  strong  heterogeneities  are  likely  to  generate  obliquely  scattered  waves  and  P/S 
conversion. 

At  the  time,  Wu  and  Flatt4  (1990)  noted  that  the  phase  screen  method  they 
used  provided  the  only  practical  means  to  numerically  simulate  the  3-D  problem.  Since 
then,  Fisk  and  McCartor  (1991)  have  developed  a  phase  screen  method  for  vector  wave 
propagation.  We  compare  simulated  TCF’s  based  on  the  method  for  vector  waves  to 
analytic  and  simulated  results  based  on  the  parabolic  approximation  to  the  scalar  wave 
equation.  We  perform  this  comparison  for  a  2-D  version  of  the  model  proposed  by 
Flatt4  and  Wu  (1988).  Although  this  study  could  be  performed  for  the  3-D  model,  for 
the  purpose  of  assessing  the  assumptions  involved  and  understanding  the  scattering 
effects,  the  2-D  model  is  adequate.  To  begin,  we  present  the  analytic  calculation.  The 
details  of  the  simulation  and  the  model  of  the  crust  and  upper  mantle  follow. 

2.1.  The  Analytic  Calculation 

Following  the  derivation  of  the  TCF’s  by  Wu  and  Flatty  (1990),  it  is  assumed 
that  the  propagation  of  the  initial  P  arrival  is  governed  by  the  scalar  wave  equation 

where  w  is  the  frequency,  ©(x)  is  the  local  P  wave  speed,  and  U  is  the  single-frequency 
P  wave.  For  weak  velocity  perturbations,  eq.  (1)  may  be  written  as 

(V*  +  Jfc*)  17  =  -2k^6n[y)U,  (2) 

where  A  =  tw/a,  a  is  the  average  wave  speed,  and  dn  =  a/a  —  1  -da/a.  Expressing 
the  wavefield  as  U  =  Uoe*,  where  Uq  is  taken  to  be  the  solution  of  the  homogeneous 
wave  equation,  one  obtains 


(V^  +  k^)  {Uork)  =  -UoiVrP  •  VrJ;  +  2kHn).  (3) 

The  solution  may  be  written  as  an  integral  equation,  and  assuming  |VV*|  is  small,  it 
takes  the  form 
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(4) 


V'(x)  =  2A;*  Uq{x)~^  d^x'  G(x,x')  Uq{x!)  5n(x'), 

where  G  is  the  Green’s  function  of  the  homogeneous  wave  equation,  and  N  is  the 
number  of  spatial  dimensions.  The  assumption  which  neglects  the  nonlinear  term 
W)  •  is  known  as  the  Rytov  approximation. 

This  expression  may  be  further  simplified  if  the  scattering  is  predominantly 
in  the  forward  direction,  t.e.  by  assuming  the  parabolic  approximation.  Let  z  be  the 
direction  of  forward  propagation  and  tt  be  the  vector  of  transverse  coordinates.  Let 
the  Green’s  function  be  expressed  in  the  form 

G(x,x')  =  I  d^-^Kr  G(Kr;  z  -  z')  (5) 

where  G{Kt',z  —  z')  is  the  Fourier  transform  w.r.t.  the  transverse  coordinates  rr  of 
G.  For  small-angle  propagation  in  N  dimensions 


G(Kr;« 


t 


2k 


Also,  assuming  a  normally  incident  plane  wave. 


(6) 


C^o(x)-^C^o(x')  (7) 

To  complete  the  analysis,  the  fiuctuations  6n  are  characterized  by  their 
power  spectrum  Pn  (K;z)  and  rms  perturbation  a  as 


(5n(x)6n(x'))  d^K  P„(K;  z)  (8) 

The  z  dependence  represents  the  possibility  that  the  power  spectrum  may  vary  with 
depth.  Now  let  ^  =  u  +  i<f>,  where  u  is  the  log  amplitude  and  ^  is  the  phase  of  the 
fiuctuating  wavefield.  Using  the  expressions  above,  the  TCF’s  for  a  2-D  medium  are 
given  by 


(u(rx)tt(O))  =  27rfc*a*  dz  j  dKr  Pn{KT,0]z)  sm^{K^z/2k)  cos{KTrT)  (9) 
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(^(rj.)^(0))  =27rA;*<7*J^  dz  j  dKj  PniKriO^z)  coB^(Kj>zf2k)  cos(iifT>’T)  (10) 


{u{rT)<f>{0))  =  nk^a^  j  dz  j  dKr  Pni^TiO'tz)  Bm{Kxz/k)  cos{KTrT)>  (11) 

In  these  expressions,  K  =  (Kr,  A*.)  and  h  is  the  depth  of  the  structure.  The  angled 
brackets  denote  an  ensemble  average  over  all  possible  realizations  of  the  medium. 
As  Wu  and  Flatt6  (1990)  have  noted,  the  wavefield  fluctuations  depend  only  on  the 
D.C.  component  of  the  heterogeneities  in  the  direction  of  propagation,  t.e.  on  the 
power  spectrum  with  K,  =  0.  Thus  the  results  are  insensitive  to  the  small-scale 
structure  in  the  vertical  direction.  Typically,  the  TCP’s  are  normalized  such  that 


(u(rr)u(0)>^  = 


{u(rr)u(0)) 

(u(0)’) 


{<f>{rT)<l>{0))N  = 


(<^(rr)(^(0)) 

{^(0)*> 


(u(rT)^(0))Ar 


(u(rr)^(O)) 


(12) 

(13) 

(14) 


These  functions  are  clearly  independent  of  the  rms  wave  speed  perturbation  a,  which 
may  be  traced  back  to  the  Rytov  approximation.  We  will  return  to  these  results 
momentarily.  First,  we  describe  how  the  simulated  results  were  obtained. 

2.2.  Phase  Screen  Simulations 

The  phase  screen  method  based  on  the  parabolic  approximation  to  the  scalar 
wave  equation  has  a  great  deal  of  similarity  to  the  calculation  just  described.  Propa¬ 
gation  between  any  consecutive  pair  of  screens  is  accomplished  by  using  the  parabolic 
free-field  Green’s  function  in  wavevector  space.  Thus  the  wavefield  at  ^'  =  z  +  Az  in 
terms  of  the  wavefield  at  z  is  given  by 


U{Kt,  z')  =  U{Kt,  z)  exp  [i{k  -  K^/2k)Az]  (15) 
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This  expression  is  now  Fourier  transformed  w.r.t.  the  transverse  component  of  the 
wavevector,  to  obtain  17 (ry, «').  The  affect  of  the  wave  speed  perturbations  are  treated 
by  multiplying  this  expression  by  the  phase  factor  exp[t^(rr)],  where 

/*+A* 

dz'  6n[TTiZ').  (16) 

Now  inverse  Fourier  transform  U{rTtz')  exp(t^(r7')]  to  obtain  (/{KrtZ'^)  on  the  other 
side  of  the  screen.  The  process  may  be  repeated  to  propagate  further.  Once  the  wave 
has  been  propagated  to  the  surface,  TCF’s  may  be  computed  numerically. 

The  fundamental  difference  between  this  approach  and  the  analytic  calcu¬ 
lation  is  the  treatment  of  diffraction  effects.  The  phase  screen  solution  divides  the 
medium  into  segments  such  that  the  geometrical  optics  regime  applies  between  con¬ 
secutive  screens.  Thus  instead  of  multiplying  the  free-held  solution  by  an  amplitude 
and  a  phase,  exp(u  + 1^),  it  is  multiplied  by  only  a  phase,  exp(t^).  Diffraction  effects 
(variations  in  the  amplitude)  are  accumulated  as  the  wave  propagates  through  many 
screens.  The  diffracted  waveheld  in  this  case  may  depend  nontrivially  on  the  strength 
of  the  perturbations.  The  Rytov  solution,  on  the  other  hand,  computes  the  diffraction 
effects  over  the  full  propagation  distance  in  a  single  step.  Unfortunately,  to  do  so,  the 
Rytov  approximation  leads  to  a  dependence  on  the  rms  perturbations  that  is  simply 
a  multiplicative  constant. 

The  phase  screen  method  for  vector  waves  has  two  additional  features.  First, 
P/S  conversion  is  treated  as  described  in  the  first  section  of  the  paper.  Second,  the 
S  and  P  waves  are  propagated  between  the  screens  using  free-held  Green’s  functions 
of  the  form  exp  [t(/c^  -  K^y/^Az],  where  k  —  wjaioi  the  P  wave  and  k  =  u;//5  for 
the  S  wave.  The  parabolic  Green’s  function  may  be  recovered  from  this  expression  by 
expanding  the  argument  of  the  exponential  to  first  order  for  Kt  <  k. 

To  generate  random  realizations  of  the  phase  screens,  we  relate  the  power 
spectrum  of  the  phase,  denoted  to  the  power  spectrum  of  the  perturbations 

7n(-Kr, Kg)  (Knepp,  1983;  Martin  and  Flatte,  1988)  via 


P^{KT)cKPn[KT,0). 


(17) 


If  the  spectra  are  normalized  such  that 

1  =  jdKr  P^{Kt)  =  f  dKrdKg  P„(Kr,if,),  (18) 
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then  the  variance  of  the  phase  is  given  by 


al  =  a^k'^Az  j  dKr  P„(ifT,0).  (19) 

Each  phase  screen  is  created  by  filtering  computer-generated  Gaussian  raiidom  num¬ 
bers  by  the  square  root  of  the  phase  spectrum.  The  rms  phase  of  each  screen  is 
normalized  by  a^. 

2.3.  The  Two-Layer  Power-Law  Model 

The  model  proposed  by  Flatt4  and  Wu  (1988)  for  the  crust  and  upper  mantle 
beneath  NORSAR  is  described  by  two  overlapping  layers  whose  heterogeneities  are 
characterized  by  power-law  spectra.  A  schematic  of  the  model  is  shown  in  Figure  2. 
The  spectra  of  the  heterogeneities  are  band  limited.  The  lower  cutoff  is  determined 
by  the  fact  that  wavenumbers  below  27r/jD  cannot  be  sampled  by  the  finite  array 
with  aperture  D.  For  NORSAR  the  array  aperture  is  D  =  110  km.  The  subarray 
beamforming  process  determines  the  upper  cutoff;  for  NORSAR  it  is  2it/5.5  km. 
The  upper  layer  is  described  by  a  flat  spectrum  down  to  a  depth  of  200  km.  This 
represents  the  fact  that  the  corner  wavenumber  of  this  spectrum  is  greater  than  the 
upper  cutoff.  Since  the  corner  wavenumber  is  inversely  related  to  the  correlation 
length  of  the  medium,  this  suggests  that  there  is  considerable  structure  at  distance 
scales  <  1  km.  The  lower  layer  ranges  from  a  depth  of  15  km  down  to  250  km.  It 
is  thought  to  consist  of  predominantly  smooth  large-scale  structure,  which  may  be 
modeled  by  an  exponential  autocorrelation  function.  The  corresponding  spectrum  in 
N  dimensions  is  given  by  a^/(l  +  where  a  is  the  correlation  length. 

Evidence  suggests  that  the  correlation  length  is  larger  than  the  lower  cutoff.  Thus  the 
band  limited  spectrum  may  be  approximated  by  a  pure  power-law  spectrum. 

The  relative  strengths  of  the  spectra  are  assumed  to  be  equal  at  if  =  0.31  km“^, 
i.e.  at  a  length  scale  of  20  km. 

3.4.  Results 

The  two  dimensional  version  of  the  model  is  obtained  by  setting  N  =  2. 
The  phase  screen  results  were  simulated  using  a  total  of  24  phase  screens;  4  were 
used  for  the  lowest  segment,  18  for  the  overlapping  segment,  and  2  in  the  uppermost 
segment.  For  the  overlapping  segment  the  phase  screens  were  generated  as  the  sum 
of  realizations  of  the  two  spectra.  This  assumes  the  two  layers  are  statistically  inde¬ 
pendent.  The  incident  wave  in  all  of  these  calculations  was  a  1  Hz  normally  incident 
plane  P  wave.  To  obtain  statistical  results  from  the  simulations,  100  realizations  were 
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averaged.  The  TCP’s  for  the  vertical  component  of  the  simulated  vector  wave  are 
compared  with  the  simulated  and  analytic  scalar  wave  results. 

Figures  3-5  show  the  TCP’s  for  this  model  with  <t  =  0.1%,  0.5%,  1.0%.  TJie 
solid,  dashed,  and  dotted  curves  represent  the  simulated  vector  wave,  simulated  scalar 
wave,  and  analytic  results,  respectively.  For  wave  speed  perturbations  <  0.5%  the 
comparison  of  the  three  calculations  is  excellent.  For  a  =  1.0%,  however,  the  simulated 
results  begin  to  differ  noticeably  from  the  analytic  results.  For  a  =  2.0%,  5.0%  (Figure 
6  and  7),  there  is  a  slight  difference  between  the  rate  at  which  the  simulated  vector 
and  scalar  TCP’s  decorrelate.  The  diflference  is  on  the  order  of  a  couple  of  kilometers. 
(Note  that  only  the  simulated  results  are  shown  in  Figures  6  and  7,  and  the  axes  are 
scaled  differently  than  in  Figures  3-5.) 

Figures  6  and  7  also  show  TCP’s  of  the  horizontal  component  of  the  scattered 
vector  wave.  It  is  not  clear  whether  the  transverse  component  TCP’s  for  a  direct  P 
arrival  can  be  put  to  significant  use;  in  practice  the  signal-to-noise  ratio  may  be  too 
small.  We  have  included  these  results  to  illustrate  that  the  phase  screen  method  for 
vector  waves  may  be  used  to  analyze  three-component  data,  perhaps  for  arrivals  other 
than  the  incident  P,  PKP  or  PKIKP  phases.  The  additional  information  contained  in 
the  transverse  components  could  possibly  lead  to  tighter  constraints  on  the  models  of 
the  crust  and  mantle  beneath  seismic  arrays  or  single  three-component  stations. 

There  are  clearly  measurable  differences  in  the  results  depending  on  the 
scattering  strength  of  the  medium.  The  simulations  verified  the  analytic  results  for 
sufficiently  weak  structures,  but  showed  that  the  wavefield  decorrelated  more  rapidly 
than  the  analytic  results  indicated  for  stronger  media.  This  is  intuitively  as  one 
would  expect;  the  wavefield  in  a  homogeneous  medium,  would  be  totally  coherent. 
The  fact  that  the  simulated  vector  wave  decorrelated  on  slightly  shorter  distances 
than  the  simulated  scalar  wave  may  be  attributed  to  the  additional  scattering  modes 
in  the  vector  case.  Only  ais  the  strength  of  the  perturbations  increased  above  2%  were 
P/S  conversion  effects  noticeable  in  the  TCP’s.  Qualitatively,  these  results  would 
also  be  observed  for  the  three-dimensional  problem.  Quantitatively,  however,  the 
generalization  is  not  readily  apparent.  Further  array  analysis  using  phase  screen 
simulations  of  elastic  waves  in  3-D  media  should  be  performed. 

3.  ACCURACY  OF  DIRECT  VERSUS  SCATTERED  PHASES 

Accurate  yield  estimates  based  on  seismic  techniques  depend  on  the  uncer¬ 
tainties  in  seismic  magnitudes.  Among  other  cispects,  the  observed  seismic  magnitudes 
are  subject  to  the  variability  of  the  intermediate  receiver-source  geology,  and  the  num¬ 
ber  and  quality  of  observing  stations.  Richards  (1988)  has  noted  that  teleseismic  body 
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waves  do  not  relate  reliably  to  the  yield  unless  the  data  is  processed  from  many  sta¬ 
tions.  Other  phases  such  as  Lg  or  P-coda  often  provide  excellent  estimates  of  yield 
for  data  obtained  at  only  one  or  two  stations,  but  depend  on  the  crustal  structure, 
propagation  distance,  etc.  In  a  study  by  Nuttli  {1986a),  22  Nevada  Test  Site  (NTS) 
nuclear  explosions  in  hard  rock  were  analyzed.  The  magnitude-yield  relation  for 
mb{Lg),  using  data  from  only  2  or  3  stations,  versus  the  announced  yield,  exhibited  a 
remarkably  small  variance.  Nuttli  (1986b)  also  analyzed  Lg  data  taken  at  NORSAR 
for  23  Semipalatinsk  explosions.  Richards  (1988)  has  noted  the  Nuttli’s  Lg  data  do 
not  exhibit  the  same  degree  of  consistency  as  that  of  published  P-coda  data  (Gupta 
et  al.,  1985),  observed  at  a  single  teleseismic  station,  for  the  same  events.  The  con¬ 
sistency  of  the  data  were  determined  by  comparing  with  an  mt(P),  estimated  from 
ISC  data  (Marshall  et  al.,  1984)  using  hundreds  of  stations.  Richards  has  suggested 
that  the  considerably  longer  propagation  distances  (>  2000  km)  may  be  responsible 
for  the  inferior  results  based  on  mi,{Lg)  for  the  Semipalatinsk  events. 

These  examples  illustrate  how  the  number  of  observing  stations,  geology  and 
propagation  distances  influence  which  portion,  or  phase,  of  the  seismogram  provides 
the  best  estimate  of  yield.  Ideally,  large  sets  of  physical  data  could  be  used  to  calibrate 
the  seismic  stations  and  determine  which  magnitudes  or  weighted  combination  of 
magnitudes  provide  the  most  reliable  yield  estimates.  Unfortunately,  relevant  seismic 
data  are  not  always  abundant,  particularly  at  new  stations.  A  practical  approach 
to  supplement  the  information  gained  from  physical  data  is  to  numerically  simulate 
seismic  data. 

As  a  first  step  towards  understanding  the  scattering  affects  responsible,  we 
have  simulated  elastic  wave  propagation  in  uniform  isotropic  random  media.  We 
generated  a  large  statistical  ensemble  of  data  (20  realizations  of  each  model  by  16 
synthetics,  at  each  25  km  interval  out  to  200  km),  and  computed  the  variances  of 
the  direct  P  wave  peak-to-peak  amplitude  and  the  rms  scattered  amplitude  of  the 
transverse  component  in  the  velocity  window  between  the  P  and  S  wave  speeds.  It 
is  the  first  such  study  of  this  type  of  which  we  are  aware.  It  is  made  possible  by  the 
tremendous  efficiency  of  the  numerical  propagation  algorithm  we  use.  The  CPU  time 
required  to  run  each  model  was  under  3  hrs.  on  an  ELXSI  6400.  To  run  this  problem 
on  a  CRAY-2  would  take  under  20  minutes. 

The  random  media  are  conveniently  characterized  by  their  autocorrelation 
functions,  which  measure  the  similarity  of  nearby  points  in  the  medium.  In  this 
study,  we  focused  on  isotropic  exponential  {N{r)  =  e“'‘/“)  and  0th  order  Von  Kar- 
man  (self-similar)  (iV(r)  =  Ko{r/a))  functions,  where  r  is  the  spatial  offset  and  a  is 
the  correlation  length  (e.g.,  Chernov,  1960;  Tatarskii,  1961).  For  all  of  the  media 
considered  here,  we  assumed  mean  P  and  S  wave  velocities  of  a  =  6.0  km/s  and 
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=  3.5  ktn/s,  respectively.  We  also  assumed  that  the  ratio  of  P  to  S  wave  speeds 
was  fixed  throughout  the  medium.  Let  the  rms  velocity  perturbation  be  denoted  by 
a  =  (^a/a)rm«  =  The  random  medium  models  considered  in  this  study 

are  listed  in  Table  1. 


Table  1.  List  of  Models. 


Model 

Corr.  Function 

a 

a  (km) 

1 

Exponential 

2% 

5.0 

2 

Self-Similar 

2% 

5.0 

3 

Exponential 

2% 

2.5 

4 

Self-Similar 

2% 

2.5 

5 

Exponential 

5% 

5.0 

6 

Self-Similar 

5% 

5.0 

7 

Exponential 

5% 

2.5 

8 

Self-Similar 

5% 

2.5 

Statistically  representative  realizations  of  the  velocities  may  be  generated  by 
the  well  established  technique  of  filtering  computer-generated  Gaussian  random  num¬ 
bers,  corresponding  to  each  point  on  a  discrete  lattice  of  the  medium,  with  the  square 
root  of  the  desired  power  spectrum.  The  power  spectrum  is  the  Fourier  transform 
of  the  autocorrelation  function.  In  two  dimensions,  the  power  spectra  corresponding 
to  the  exponential  and  self-similar  functions  are  a*(l  +  and  o*(l  + 

respectively.  The  media  were  normalized  by  their  standard  deviation  a.  For  this 
study,  the  number  of  gridpoints  used  to  sample  the  medium  was  512  by  2048,  and  the 
sampling  interval  was  0.1  km  in  each  direction. 

Once  the  velocities  for  a  particular  medium  have  been  generated,  the  phase 
factors  needed  in  the  phase  screen  algorithm  were  obtained  as  follows.  The  medium 
is  divided  into  intervals,  by  planes  of  constant  z  (referred  to  as  screens  in  much  of  the 
literature),  such  that  the  awicumulated  phase  between  the  screens  may  be  computed 
from  geometrical  optics  principles.  The  phase  for  the  P  wave,  computed  by  integrating 
(or  for  the  discrete  case,  summing)  the  P  wave  speed  perturbations  along  the  direction 
of  propagation,  is  given  by  eq.  (16).  There  is  an  analogous  expression  for  the  S  wave 
phase  factor.  In  practice,  40  phase  screens  were  used  in  each  run. 
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A  Ricker  wavelet  was  used  as  the  initial  pulse  in  our  study.  The  time  source 
function  of  a  Ricker  wavelet  and  its  Fourier  transform  are  given  by 


F(t)  =  (1  -  (27r/o(t  -  to))V2)  exp  (-(27r/o(t  -  to))*/^)  .  (20) 

27ritof)  ,  (21) 

where  fo  is  the  peak  frequency  and  to  defines  the  origin  of  time.  The  value  /o  =  1 
Hz  vfdis  used  for  all  of  the  simulations.  The  maximum  frequency  sampled  was  5 
Hz,  where  the  norm  of  the  spectrum  has  dropped  off  to  ~  10“^°.  The  number  of 
positive  frequencies  sampled  was  256;  the  remaining  frequency  data  were  zero  padded 
to  minimize  wrap-around  effects  when  numerically  Fourier  transforming  to  the  time 
domain. 


Sixteen  synthetic  seismograms  were  generated  at  25  km  intervals  out  to  200 
km  for  twenty  realizations  of  each  random  medium  model  considered.  The  synthet¬ 
ics  were  evenly  spaced  in  the  horizontal  direction.  Figures  8-15  show  characteristic 
synthetics  at  50  km  intervals  for  the  eight  models.  Using  a  total  sample  of  320  seis¬ 
mograms  at  each  interval,  we  computed  the  first  and  second  statistical  moments  of 
the  vertical  peak-to-peak  P  wave  amplitude  (mt(P)),  and  the  rms  amplitude  of  the 
transverse  component,  in  the  time  window  between  2:/a:  -f  l//o  to  zfp.  The  time 
window  corresponds  more  closely  to  that  of  P-coda,  however,  we  have  denoted  this 
rms  amplitude  by  mi,{Lg)  since  it  is  part  of  the  transverse  component.  (We  hope  that 
the  serious  seismologist  will  not  hold  this  notation  against  us.  A  layered  structure, 
necessary  to  actually  produce  Lg,  is  not  present  in  our  models.  We  have  also  misused 
the  definition  of  mj,  which  is  customarily  used  to  denote  the  logarithm  of  an  ampli¬ 
tude.  These  expressions  should  be  thought  of  only  as  amplitudes  of  direct  (vertical) 
and  scattered  (horizontal)  phases.) 

In  Figures  16-23  we  have  plotted  the  standard  deviations,  given  as  percent¬ 
ages  relative  to  the  means.  These  plots  show  that  near  the  source  the  direct  P  wave 
amplitude  has  far  less  variation  than  the  scattered  phase,  which  is  only  beginning  to 
form.  However,  at  some  distance  away,  between  50  km  to  150  km  from  the  source 
for  these  models,  the  scattered  phase  has  less  variation.  In  all  cases,  the  cross-over 
distance  was  greater  for  the  self-similar  media  than  the  exponential  media  with  the 
same  parameters.  This  result  may  be  attributed  primarily  to  the  rate  at  which  the 
P  wave  standard  deviation  changes  with  distance,  and  from  one  model  to  next.  In 
fact,  perhaps  the  most  interesting  result  of  this  study  is  that  the  standard  deviation 
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of  the  scattered  phase  is  so  insensitive  to  the  model.  At  200  km  the  relative  standard 
deviation  of  the  scattered  phase  is  between  0.11  and  0.12  for  six  of  the  eight  models. 
The  other  two  values  are  0.13  and  0.14  for  models  1  and  5,  respectively.  In  contrast, 
at  200  km  the  relative  standard  deviations  of  the  P  wave  range  from  0.14  for  model  4 
to  0.41  for  model  1. 

To  interpret  these  results,  scattering  theory  provides  some  insight.  Wu 
(1990)  has  pointed  out  that  large-scale  heterogeneities  have  the  greatest  affect  on  the 
distribution  of  forward  scattered  waves.  Thus  it  is  not  surprising  that  the  media  with 
the  greatest  power  at  low  wavenumber  (or  equivalently,  on  large  distance  scales)  should 
produce  the  largest  variations  in  the  direct  P  wave.  Models  with  5%  perturbations 
lead  to  greater  variations  in  the  P  wave  than  those  with  2%.  Also,  since  the  self¬ 
similar  media  have  more  small-scale  structure  than  the  exponential,  but  both  are 
normalized  to  have  the  same  total  power  a,  the  power  at  low  wavenumber  is  greater 
for  the  exponential  medium.  Our  results  indicate  that  exponential  media  lead  to 
larger  variations  in  the  P  wave  amplitude  than  do  the  corresponding  self-similar  media, 
keeping  everything  else  fixed. 

Media  with  larger  correlation  lengths,  for  all  other  parameters  fixed,  lead 
to  larger  variations  in  the  P  wave  as  well.  One  might  argue  that  this  result  could  be 
can.^^d  in  part  by  the  fact  that  we  have  averaged  over  twice  as  many  synthetics  per 
correlation  length  when  a  =  2.5  km  as  when  a  =  5.0  km.  To  explore  this  possibility,  we 
have  compared  the  results  when  only  half  as  many  samples  per  correlation  length  were 
used  for  the  a  =  2.5  km  cases.  These  results  are  represented  by  the  open  symbols  in 
Figures  18,19,22,23.  The  difference  in  the  relative  standard  deviations  are  negligible. 

Wu  (1990)  has  also  noted  that  the  small-scale  heterogeneities  are  responsible 
for  scattering  at  large  angles,  and  hence  the  generation  of  coda.  Although  the  relative 
standard  deviations  of  the  scattered  phase  do  not  vary  much  from  one  model  to  the 
next,  they  are  in  general  the  smallest  for  the  self-similar  media.  The  presence  of 
considerable  small-scale  heterogeneities  produces  a  more  uniform  distribution  in  the 
scattered  waves. 

Simulations  for  more  realistic  layered  earth  models  are  needed  before  further 
conclusions  may  be  drawn  concerning  the  variability  of  body  waves,  surface  waves, 
regional  phases  and  P-coda.  Our  study  shows  that  useful  information  may  be  gained 
from  such  simulations.  It  will  be  necessary  to  use  an  efficient  propagation  algorithm, 
such  as  the  phase  screen  method  used  here,  in  order  to  generate  a  sufficiently  large 
synthetic  database  needed  for  this  type  of  study. 
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Figure  1.  The  multiple  phase-icreen  method  replaces  (a)  each  segment 
of  the  heterogeneous  medium  with  (b)  a  uniform  segment  and 
a  phase  screen.  The  accumulated  position-dependent  phase  is 
projected  onto  screen  3. 


Figure  2.  Schematic  diagram  of  .Flatt4  and  Wu’i  twa>layer  power- 
law  medium  for  modeling  the  heterogeneities  beneath  NOR- 
SAR.  The  power  spectra  for  the  different  layers  are  shown 
on  the  right  (from  Flattd  and  Wu,  1098). 
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Figure  3.  IVansverse  coherence  functiong  of  the  log  amplitude,  phase 
and  cross  correlation  computed  for  a  P  wave  transmitted 
through  the  two-layer  power-law  model  for  an  rms  velocity 
perturbation  of  0.1%.  The  solid,  dashed,  and  dotted  curves 
represent  the  simulated  vector  wave,  simulated  scalar  wave, 
and  analytic  scalar  wave  results. 
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Same  as  Figure  3,  but  for  an  rms  velocity  perturbation 
0.5%. 
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Figure  5. 


Same  as  Figure  but  for  an  rms  velocity  perturbation  of 

1.0%. 
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Figure  6.  TCF*s  for  the  vertical  component  of  the  simulated  vector 
wave  (solid  curve)  and  the  simulated  scalar  wave  (dashed 
curve)  are  on  the  left.  On  the  right  are  TCF’s  computed 
for  the  transverse  component  of  the  simulated  vector  wave. 
The  rms  velocity  perturbation  was  2.0%. 
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Figure  7.  Same  as  Figure  6,  but  for  an  rms  velocity  perturbation  of 
5.0%. 
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trace)  components  are  plotted  in  pairs. 
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Figure  10. 
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Synthetic  seismograms  for  model  3  at  60  km  intervals  out 
to  200  km.  The  vertical  (upper  trace)  and  horizontal  (lower 
trace)  components  are  plotted  in  pairs. 
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Figure  12.  Synthetic  geismograms  for  model  5  at  50  km  intervals  out 
to  200  km.  The  vertical  (upper  trace)  and  horizontal  (lower 
trace)  components  are  plotted  in  pairs. 
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Figure  13.  Synthetic  seismograms  for  model  6  at  50  km  intervals  out 
to  200  km*  The  vertical  (upper  trace)  and  horizontal  (lower 
trace)  components  are  plotted  in  pairs. 
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Figure  14.  Synthetic  seitmogrami  for  model  7  at  60  km  intervali  out 
to  200  km.  The  vertical  (upper  trace)  and  horizontal  (lower 
trace)  component!  are  plotted  in  pain. 
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Synthetic  seismograms  for  model  8  at  50  km  intervals  out 
to  200  km.  The  vertical  (upper  trace)  and  horizontal  (lower 
trace)  components  are  plotted  in  pairs. 
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Figure  16.  Model  1  standard  deviations,  given  as 
peak-to-peak  P  wave  amplitudes 
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Figure  17.  Model  2  standard  deviations,  given  as  percentages  relative  to  the  means,  of  the  direct 
peak-to-peak  P  wave  amplitudes  <Tm*(p),  aud  the  rms  transverse  amplitudes 
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Figure  21.  Model  6  standard  deviations,  given  as  percentages  relative  to  the  means,  of  the  direct 
peak-to-peak  P  wave  amplitudes  and  the  rms  transverse  amplitudes 


REVIEW  OF  ATTENUATION  IN  SALT  AT  MODERATE  STRAINS 


W.  R.  Wortman  and  G.  D.  McCartor  ^ 


Abstract.  Experiments  which  reflect  the  attenuation  of  propagating  pulses  in  salt  in  the 
moderate  strain  regime  of  10’^  to  10"^,  which  corresponds  roughly  to  ranges  of  100  to  10,000 
meters  from  an  explosion  with  a  yield  of  1-kt,  are  reviewed.  A  transition  from  nonlinear  to  linear 
behavior  occurs  in  this  interval.  This  regime  is  important  for  monitoring  of  nuclear  test  treaties 
since  models  for  linear  source  functions  are  normally  defined  inside  it.  Salt  is  of  interest  since  it 
can  be  readily  mined  to  produce  decoupling  cavities.  Data  from  explosive  sources,  n  tclear  and 
chemical  field  and  small  scale  laboratory  tests,  resonant  bars  and  ultrasonic  pulse  me  . hods  are 
summarized.  The  experiments  are  diverse  in  their  character  and  frequency  and  no  single 
experiment  covers  the  nonlinear-linear  transition.  However,  the  totality  suggests  that  attenuation 
in  salt  does  decrease  dramatically  over  the  moderate  strain  regime.  A  full  physical  description 
does  not  exist  although  shear  failure  or  yielding  can  account  for  some  effects. 

Introduction 

Seismic  monitoring  of  underground  nuclear  tests  requires  that  properties  of  distant  seismic 
signals  be  related  to  the  yield  of  the  explosive  source.  In  general,  providing  this  relation  requires 
finding  the  material  response  to  explosive  loading  under  an  extreme  range  of  conditions.  The 
rock  next  to  the  explosive  device  is  vaporized  out  to  about  10  meters  for  a  1-kt  explosion. 
Somewhat  beyond  this  the  rock  is  melted  and  crushed.  Out  to  a  range  of  many  tens  of  meters  the 
rock  is  visibly  cracked.  All  these  effects  are  highly  nonlinear.  Fortunately,  it  is  not  necessary  to 
provide  first  principles  calculations  of  this  behavior  in  order  to  monitor  nuclear  tests.  Subsurface 
ground  data  have  been  taken  which  can  be  used  to  define  the  pulse  from  explosions  in  a  variety 
of  media  outside  the  highly  nonlinear  regime.  These  data,  along  with  some  supporting 
theoretical  calculations  for  interpolation,  can  be  used  to  define  an  equivalent  seismic  source  at  the 
range  where  the  data  were  taken.  This  bypasses  many  uncertainties  from  complex  material 
behavior  under  extremes  of  temperature  and  pressure.  It  is  common  to  assume  that  the  resulting 
source,  which  is  typically  defined  at  a  few  hundred  meters,  can  be  propagated  further  to  seismic 
receivers  assuming  that  the  medium  is  strictly  linear,  although  not  elastic  with  attenuation 
described  by  a  linear  Q  operator,  while  taking  into  account  the  intervening  geologic  structure. 
With  this,  the  computed  waveform  properties  can  be  compared  with  those  observed  to  estimate 
yield  and  depth.  This  technique  is  the  fundamental  approach  to  monitoring  nuclear  test  treaties. 
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As  we  shall  see  in  this  article,  il^ere  is  strong  evidence  of  nonlinear  attenuation  in  the  ranee  of 
moderate  strains  (sometimes  called  high-strains),  which  we  define  here  as  strains  from  10"^  to 
10"^  which  occurs  at  roughly  100  to  10,000  m  from  a  1-kt  explosion,  which  extends  well 
outside  the  usual  region,  or  "elastic  radius,"  where  equivalent  seism?';  sources  are  defined.  The 
effect  of  this  nonlinear  attenuation,  for  which  the  attenuation  depem's  on  the  amplitude  of  the 
disturbance  and  not  just  the  material  properties,  between  ranges  showing  highly  nonlinear 
structural  changes  and  ranges  with  linear  strains  levels  (roughly  10"^)  is  not  well  understood. 
There  is  a  substantial  body  of  data  on  attenuation  in  salt  which  can  be  brought  to  bear  on  this 
question.  Salt  is  also  particularly  interesting  because  it  is  relatively  easily  mined  to  produce 
cavities  which  could  serve  to  decouple  clandestine  tests  [Evemden  et  al..  1 986].  For  both  of 
these  reasons,  salt  is  a  medium  for  which  the  seismic  response  is  important.  In  this  article  we 
shall  indicate  the  nature  and  results  of  experiments  which  have  been  done  in  salt.  There  are  three 
types  of  experiments  which  have  provided  useful  attenuation  data  in  salt.  These  are:  explosive 
source  generated  pulse  propagation,  either  in  the  field  or  the  laboratory';  ultrasonic  pulse 
propagation;  and  resonant  bar  excitation. 

The  explosive  source  technique  is  obviously  most  directly  related  to  the  question  at  hand.  It 
consists  of  placing  a  set  of  motion  sensors  at  various  ranges  and  orientations  relative  to  a  source 
to  cover  the  amplitudes  of  interest.  As  a  practical  matter,  there  are  limits  on  the  number  of 
sensors  which  can  be  used  due  to  the  available  volume  of  uniform  medium.  The  result  is  a  set  of 
time  domain  pulses  measured  at  a  discrete  set  of  sensors  which  can  be  ordered  according  to  range 
from  the  source  to  the  sensor.  When  the  effects  of  orientation,  instrument  response  and  spherical 
divergence  are  removed,  the  change  of  the  pulse  with  range  reflects  the  results  of  any  attenuation. 
For  linear  attenuation  described  by  a  Q  operator,  the  multiplicative  change  in  amplitude  in  the 
frequency  domain  with  a  range  change,  5r,  goes  like 

exp(-(j)5r/2cQ)  (1) 

where  c  is  the  phase  speed.  For  linear  mechanisms,  Q*^  can  be  interpreted  as  the  fractional 
energy  loss  per  cycle  over  47i  or  the  tangent  of  the  phase  angle  between  the  stress  and  strain. 
Through  the  above  relation,  generally  referred  to  as  using  the  spectral  ratio  method,  Q  is  a 
common  measure  of  attenuation  even  if  the  mechanism  is  not  known  to  be  linear.  Using  this 
measure,  a  Q  as  a  function  of  frequency  can  be  estimated  for  every  pair  of  records. 

Ultrasonic  pulse  propagation  using  a  planar  geometry  can  also  be  used  to  estimate  a  Q  I'rom 
(1).  A  sample,  usually  with  a  linear  dimension  of  a  few  centimeters  but  large  compared  with  the 
pulse  wavelength,  is  placed  between  a  pair  of  plates.  One  plate  is  pulsed  with  a  finite  source 
waveform  while  a  transducer  at  the  receiver  plate  detects  the  transmitted  pulse.  The  timing  of  the 
received  pulse  allows  determination  of  propagation  speed  while  the  spectral  ratio  of  the 
transmitted  and  received  pulses  fixes  Q.  Both  compressional  and  shear  pulses  can  be  measured 
by  this  method.  The  requirement  of  a  reasonable  apparatus  and  sample  size  limits  this  technique 
to  ultrasonic  frequencies. 
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The  resonant  or  oscillating  bar  method  does  not  involve  pulse  propagation.  It  measures  the 
response  of  a  damped  resonant  system  to  a  harmonic  driver,  for  which  the  damping  element  is 
the  material  under  study.  Typicdly  a  cylindrical  rod  approximately  10  cm  in  length  is  loaded 
with  appropriate  weights  to  give  a  resonant  frequency  as  low  as  practical,  often  a  few  hundred 
Hertz.  The  system  is  then  driven  in  torsion  or  flexure  modes  and  its  resonant  response 
determined  as  the  driver  frequency  is  scanned  across  the  spectrum.  The  resulting  response  given 
as  an  amplitude  as  a  function  of  frequency  determines  the  damping  of  the  system.  The  half 
power  bandwidth  of  the  resonance  curve  divided  by  the  resonant  frequency  is  Q‘^. 

It  must  be  pointed  out  that  the  concept  of  Q  is  confined  to  linear  attenuation.  If  there  are 
nonlinear  effects,  Q  may  well  not  be  a  robust  description  which  allows  a  meaningful  comparison 
of  different  experimental  results  or  a  useful  means  of  applying  the  results  of  experiments.  The 
danger  of  blind  use  of  experimental  Q  vill  be  illustrated  in  this  article. 

Experimental  Attenuation  Data  Available  For  Salt 

Explosive  Sources-  Nuclear 

Salmon.  The  nuclear  explosion  Salmon  (5.3-kt)  event  took  place  in  a  natural  salt  dome  in 
Mississippi  in  1964  at  a  depth  of  828  meters.  A  comprehensive  set  of  approximately  83  near¬ 
field  measurements  at  16  instrument  stations  was  planned,  both  at  surface  and  subsurface 
locations  at  ground  ranges  from  166  to  744  meters  [Ferret,  1967].  A  parallel  set  of 
measurements  was  carried  out  by  SRI  [Eisler  and  Hoffman,  1966].  A  rather  consistent  set  of 
near-field  subsurface  measurements  including  scaled  ranges  from  95  to  425  resulted. 

The  single  component  acceleration  and  velocity  instruments  provided  measurements  of  a  range  of 
peak  strains  from  about  4xl0'3  to  2x10'^  with  dominant  frequencies  from  1  to  over  10  Hz.  A 
representative  set  of  particle  velocity  records  [McCartor  and  Wortman,  1985]  is  shown  in 
Figure  1. 

The  Salmon  data  have  been  discussed  by  several  authors.  The  original  experimenter.  Ferret 
[1967],  gives  an  extensive  description  of  the  experimental  details.  He  indicates  that  the  peak 
acceleration,  velocity  and  displacement  data  all  tend  to  fall  off  in  a  fairly  uniform  manner  with 
range.  The  peak  radial  velocity  falls  off  like  range  to  the  -1.87610.049  power  as  shown  in 
Figure  2.  This  decay  rate,  which  is  well  in  excess  of  the  -1  power  for  geometrical  far- field 
elastic  fall  off,  is  an  indication  that  there  is  some  attenuation  mechanism  in  operation  but  it  is  not 
necessarily  an  indication  of  nonlinearity. 

Rogers  [1966]  described  the  Salmon  free-field  data  and  compared  the  results  with  finite 
difference  calculations.  These  comparisons  generally  show  agreement  with  the  magnitude  of  the 
pulses  but  the  waveforms  are  noticeably  different.  This  indicates  that  the  models  of  highly 
nonlinear  material  behavior  based  on  prior  experiments  in  salt  were  within  reason,  but  these 
mechanisms  generally  cease  to  be  important  at  the  ranges  of  the  Salmon  data.  This  article  was  a 
part  of  a  Journal  of  Geophysical  Research  issue  [Volume  71,  No.  14, 1966]  with  a  sequence  of 
articles  on  the  Salmon  event. 
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McCartor  and  Wortman  [1985]  analyzed  a  selected  sequence  of  the  Salmon  free-field  records 
for  the  purpose  of  finding  clear  evidence  of  nonlinearity  with  this  data  set  alone.  It  was  found 
that  these  data  require  an  effective  Q  which  increases  significantly  with  increasing  frequency  and 
which  has  a  value  of  about  5  in  the  dominant  frequency  range  from  5  to  10  Hz.  This  is 
consistent  with  Trulio  [quoted  by  Larson,  1982]  indicating  that  the  Salmon  data  for  decay  of 
peak  velocity  with  range  are  consistent  with  an  effective  Q  of  about  3  in  the  0.5  to  5  Hz  range  for 
strains  near  10*^.  Figure  3  shows  a  typical  effective  Q  found  using  the  spectral  ratio  of  free-field 
records  from  166  and  660  meters  from  the  shot  point.  Examination  of  a  series  of  Q  estimates 
from  5  record  pairs  over  this  same  range,  shown  in  Figure  4,  does  not  show  clear  evidence  of 
decreasing  attenuation  with  range.  The  extreme  range  record  pairs  attenuation  estimates  appear  to 
differ  by  about  20%,  decreasing  with  range,  but  this  is  the  same  level  as  the  uncertainties.  The 
Salmon  elastic  precursor  which  leads  the  Salmon  pulses  has  been  accounted  for  by  McCartor  and 
Wortman  [1990]  by  a  partial  shear  failure  which  permanently  reduces  the  shear  modulus  when 
the  strain  exceeds  a  threshold  of  about  10*^.  Rimer  and  Cherry  [1982]  have  shown  that  it  is 
possible  to  reproduce  much  of  the  Salmon  data,  including  the  precursor,  using  a  shear  strength 
limit  which  is  variable. 


Gupta  and  McLaughlin  [1989]  analyzed  Salmon  and  Sterling  data  and  coiicluded  that  the 
effective  Q  at  Salmon  strains  appears  to  be  mildly  strain  dependent  and  strongly  frequency 
dependent.  The  mean  apparent  Q  in  the  1-25  Hz  range  is  about  7  and  it  appears  to  increase 
mildly  with  increasing  rang^  The  result  for  attenuation  is  reflected  in  a  m^ified  Q  function 
called  Q.  This  is  defined  as  Q=  Q+f  dQ/df  which  is  a  measure  of  the  spectral  slope  change.  The 
behavior  found  is  shown  in  Figure  5.  Gupta  and  McLaughlin  argue  that  the  decrease  in 
attenuation  with  range  is  significant  and  this  indicates  that  the  behavior  is  nonlinear.  The 
attenuation  decreases  sharply  with  increasing  frequency.  Denny  [1990]  reports  that  the  source 
spectra  characteristics  of  Salmon  and  Sterling  indicate  that  the  Salmon  pulses  are  nonlinear  to 
beyond  700  m. 

Sterling.  The  nuclear  explosion  Sterling  (0.38-kt)  event  took  place  in  1966  in  the  Salmon 
cavity,  which  was  approximately  spherical  with  a  radius  of  17  meters.  This  was  the  second  half 
of  the  decoupling  test  and  the  same  instrumentation  was  used.  The  waveforms  observed  are 
generally  noisy  and  less  cohesive  than  for  Salmon.  Sterling  data  [Sisemore  et  al.,  1969]  are  at 
lower  strains  than  Salmon  due  both  to  lower  yield  and  decoupling.  Sterling  peak  velocities 
indicate  a  strain  range  of  3x10’^  to  7x10'^.  Some  analysis  of  these  data  by  Springer  et  al. 
[1968]  suggested  that  there  was  significant  near-field  attenuation.  A  more  recent  analysis  by 
Glenn  et  al.  [1987]  corrects  some  errors  in  the  previous  work  and  indicates  that  the  observed 
Sterling  near-field  pulses  are  in  good  agreement  with  elastic  theory.  Gupta  and  McLaughlin 
[1989]  used  a  spectral  ratio  method  to  determine  the  attenuation  over  sensor  pairs.  They  find  that 
the  average  Q  is  approximately  200  to  400  and  shows  no  evident  dependence  on  frequency  or 
range.  These  last  two  papers  suggest  that  the  Sterling  strains  have  reached  a  transition  to  a  linear 
low  attenuation  at  small  strains,  near  10'^.  Langston  [1983]  has  noted  that  SV  waves  generated 


by  Sterling,  apparently  due  to  induced  normal  faulting  rather  than  cavity  asymmetry,  showed 
attenuation  with  range  which  was  significantly  different  from  1/r.  He  finds  that  Qp  of  35  is 
indicated  at  strains  of  about  10'^  and  there  is  a  mild  tendency  for  Qp  to  increase  with  range,  and 
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so  with  decreasing  strain.  If  it  is  assumed  that  there  is  no  contribution  to  attenuation  from 
compression,  the  corresponding  P-wave  Q  will  be  approximately  70,  although  there  is  no  direct 
evidence  to  tie  the  two  modes  together  for  this  particular  case. 

Gnome-  The  Gnome  explosion  (3  kt)  took  place  in  the  Salado  salt  layer  at  a  depth  of  360 
meters  [Weart,  1962],  Above  the  salt  there  were  strata  to  a  depth  of  about  200  meters  consisting 
of  diverse  lithology.  Analysis  of  the  data  from  eight  sensors  sites  located  horizontally  from  the 
shot  point  at  ranges  of  62  to  477  meters  indicates  that  there  are  three  distinct  arrival  times.  The 
high  speed  initial  signal  appears  to  be  from  refraction  in  a  polyhalite  layer  below  the  shot,  the 
normal  moderate  speed  signal  then  appears  followed  by  a  plastic  wave  with  a  speed  which 
diminishes  with  decreasing  strain.  Some  of  the  acceleration  records  are  clipped  so  it  is  difficult 
to  use  the  details  of  the  waveform  to  characterize  the  attenuation.  Weart  indicates  that  the  peak 
velocity  falls  off  like  range  to  the  -3.56  power  out  to  100  meters  and  then  like  the  power  -1.36 
out  to  the  last  sensor  at  477  meters.  It  is  suggested  that  the  elastic  zone  has  been  reached  at  this 
extreme  range  but  no  evidence  beyond  the  amplitude  variation  with  range  is  given. 

Explosive  Sources-Chemical 

Gowbpy-  The  Cowboy  series  of  chemical  explosions  took  place  in  a  salt  dome  in  Louisiana 
in  1959-60  [Murphey,  1961].  The  explosions  had  a  range  of  yields  from  10  to  2000  pounds  of 
TNT,  some  of  which  were  carried  out  in  cavities  for  decoupling  tests.  Each  shot  was  sampled 
by  from  two  to  seven  particle  velocity  sensors.  The  scaled  ranges  for  the  coupled  or  tamped 
experiments  were  from  about  200  to  3000  m/kt^/^  and  the  corresponding  peak  strains  were  from 
a  few  times  10“^  to  about  10'5.  The  dominant  frequencies  were  10  to  100  Hertz.  Murphey 
observed  that  the  peak  vel^ity  with  scaled  range  (range/yield  data  all  lie  near  a  smooth  curve 
and  they  fall  off  like  r'l*o5^  as  shown  in  Figure  6.  This  indicates  first  that  the  data  scale  and 
second  that  the  material  behavior  is  not  elastic. 

Minster  and  Day  [1986]  examined  the  Cowboy  tamped  data  and  investigated  the  attenuation 
required  to  reproduce  the  observed  variation  of  peak  velocity  and  displacement  with  scaled  range. 
They  conclude  that  an  effective  Q  which  is  strain  dependent  with  the  frequency  independent  form 

1/Q  =  1/Qo  +  ye  (2) 

can  satisfactorily  reproduce  the  both  peak  velocity  and  displacement.  Here  Qq  =  100  is  the 
small  strain  Q,  e  is  the  peak  strain  and  y~  3x10^  is  an  empirical  constant,  this  form  was 
chosen  to  be  consistent  with  some  theoretical  nonlinear  mechanisms  at  high  strains  while 
reducing  to  modest  anelastic  attenuation  at  small  strains.  The  resulting  attenuation  is  consistent 
with  scaling  since  it  depends  only  on  strain  which  scales.  Note  that  ihe  effective  Q  at  the  largest 
strains  from  Cowboy  is  then  less  than  10. 

Trulio  [private  communication]  has  examined  some  scaled  Cowboy  data  and  concluded  that 
the  attenuation  is  frequency  dependent,  decreasing  with  increasing  frequency.  Q  values  of  12.5 
and  32  are  found  at  Salmon  equivalent  frequencies  of  10^/^  and  10  Hz,  found  by  scaling  the 
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Cowboy  data.  The  frequency  dependence  is  roughly  Q  oe  i/f.  He  observes  that  this  dependence 
is  inconsistent  with  linearity  and  scaling  indicating  that  Cowboy  attenuation  is  nonlinear. 

Wortman  and  McCartor  [1989]  have  used  tamped  Cowboy  records  to  attempt  to  determine  the 
character  of  the  attenuation.  They  chose  record  pairs  from  the  same  shot  and  applied  the  spectral 
ratio  method  to  determine  the  effective  Q  between  the  stations.  By  selecting  pairs  with  internal 
consistency  and  by  correcting  the  records  for  obvious  anomalies,  they  found  a  set  of  Q  estimates 
as  a  function  of  peak  strain  or  range.  The  result  is  shown  in  Figure  7  as  compared  with  the  result 
of  Minster  and  Day  [1986].  Not  surprisingly,  the  results  are  similar  and  they  show  that  across 
the  strains  available  in  the  Cowboy  data  there  is  a  substantial  decrease  in  attenuation  with 
decreasing  strain. 

Cowboy  Trails.  Cowboy  Trails  [Workman  and  Trulio,  1985]  was  a  series  of  field  tests  in  a 
salt  dome  using  chemical  explosives  of  approximately  200  pounds.  The  single  and  dual 
explosions  were  monitored  with  free-field  ground  velocity  sensors.  The  ranges  of  the  sensors, 
which  were  varied  over  the  ten  events,  covered  from  0.388  to  over  1 1.3  km/kt^/^.  They  were 
set  to  attempt  to  define  the  transition  to  clearly  linear  behavior.  Due  to  problems  with  sensors 
and  noise,  the  Cowboy  Trails  experiments  did  not  fully  accomplish  their  goals.  However,  the 
results  indicate  that  the  propagation  of  explosively  driven  pulses  in  dome  salt  remain  inelastic 
(though  not  necessarily  nonlinear)  out  to  scaled  ranges  of  1 1 .3  km/kt^/^.  The  exponent  of  peak 
velocity  power  law  decay  with  range  is  -1.46±.05  which  is  somewhat  slower  than  that  seen  in 
the  Cowboy  data  (which  were  for  higher  strains).  An  analysis  by  Trulio  [private 
communication]  suggests  that  the  scaled  Cowboy  Trails  data  show  a  decreased  attenuation  with 
range  and  frequency  which  is  not  inconsistent  with  scaled  Salmon  and  Cowboy  results.  The 
scatter  of  the  data  is  fairly  large  so  it  is  difficult  to  draw  stronger  conclusions. 

Livermore  1982  Small  Scale  Experiments.  The  experiments  of  Larson  [1982]  for  small 
chemical  explosions,  yields  of  0.63  to  291  kJ,  in  pressed  salt  have  provided  pulses  over  scaled 
ranges  from  approximately  10  to  250  m/kt^/^.  The  dominant  range  of  frequencies  covered  was 
from  about  10^  to  10^  Hz  and  the  ratio  of  peak  particle  velocities  to  compressional  sound  speed 
(which  is  comparable  with  the  strain)  went  from  about  10*1  to  less  than  lO'^.  Data  from  three 
sensors  at  increasing  ranges  for  a  single  shot,  taken  in  pairs,  indicate  increasing  values  of  Q 
from  12  to  25  with  increasing  range  for  ranges  from  30  to  70  m/ktl/3.  This  would  suggest  that 
the  response  was  nonlinear.  Another  experiment  consisting  of  a  simultaneous  pairs  of  shots, 
was  used  as  a  direct  superposition  experiment.  It  was  found  that  the  resulting  response  was 
consistent  with  that  determined  by  linear  addition  of  the  two  pulses  at  a  range  of  168  m/ktl/3  as 
would  be  expected  from  a  linear  medium.  Still  it  is  not  clear  just  how  nonlinear  effects  would  be 
manifest  in  this  experiment  of  rather  narrow  pulses  without  knowing  the  character  of  any 
nonlinear  behavior.  That  is,  the  apparent  agreement  with  superposition  for  pulses  with  large 
strains  may  very  well  not  directly  negate  the  possibility  of  any  sort  of  nonlinear  behavior.  Still, 
on  the  basis  of  this  experiment,  Larson  concludes  that  the  propagation  is  linear  beyond  168 
m/ktl/3.  He  also  indicates  that  the  attenuation  is  strongly  inelastic  inside  this  range  and  that  the 
magnitude  of  this  inelasticity  decreases  with  both  range  and  frequency.  By  variation  of  the 
confining  pressure,  he  finds  that  the  propagation  is  independent  of  confinement  at  least  up  to 
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32  MPa.  Larson  combined  his  data  with  that  from  other  salt  experiments  to  further  extend 
Trulio's  [1978]  assertion  that  peak  velocity  data  for  explosively  driven  pulses  in  salt  scales 
remarkably  well.  This  is  shown  in  Figure  8  which  contains  data  from  experiments  with  over  10 
orders  of  magnitude  in  yield.  Larson  also  provides  a  direct  comparison  of  scaled  waveforms  for 
his  experiment  and  a  Salmon  pulse  which  shows  a  significant  similarity. 

Livermore  1987  Small  Scale  Experiment.  In  1987  Larson  [1989]  carried  out  an  additional 
laboratory  experiment  using  a  0.622  kJ  chemical  explosion  in  mined  dome  salt.  The  cylindrical 
sample  was  sliced  into  layers  perpendicular  to  the  axis.  Sixteen  sensors  were  then  inserted  at 
eight  ranges  and  the  layers  were  cemented  together.  A  explosive  source  tamped  in  a  hole  on  one 
side  then  produced  peak  particle  velocities  from  1  m/s  to  4  xlO'^  m/s.  The  experiment  was  not 
very  successful.  A  variety  of  problems  with  the  samples  caused  irregularities  in  the  waveforms 
which  made  analysis  difficult.  Wortman  and  McCartor  [1989,  Appendix  C]  have  reproduced  the 
waveforms  and  attempted  to  find  the  attenuation.  Aside  from  an  apparent  increase  of  Q  with 
frequency,  little  can  be  concluded. 

Nonexplosive  Sources 

Ultrasonic  Pulse  Attenuation.  New  England  Research  (NER)  laboratory  ultrasonic  pulse 
propagation  experiments  [Coyner,  1987]  used  strains  from  less  than  10‘6  to  more  than  10*5. 
Compressional  and  shear  ultrasonic  pulses  consisting  of  about  two  cycles  at  100-200  kHz  were 
propagated  through  samples.  Attenuations  were  calculated  using  a  spectral  ratio  technique  with 
an  aluminum  sample  used  for  calibration.  Variation  of  the  attenuation  with  peak  strain  amplitude 
and  confining  pressure  were  determined.  For  dome  salt  it  was  found  that  over  a  strain  range  of 
5x10"^  to  3x10’5  and  for  confining  loads  of  either  0.1  or  1  MPa,  the  P-wave  attenuation  is 
nearly  constant  and  can  be  described  by  a  Q  of  about  20.  For  S  waves,  the  Qn  is  also  nearly 
constant  and  it  has  a  value  of  about  60.  The  results  are  shown  in  Figure  9.  mth  the  possible 
exception  of  the  largest  strains  for  S  waves,  there  is  no  particular  evidence  of  nonlinearity  in 
these  data  alone  although  the  attenuation  is  large.  The  P-wave  attenuation  is  about  a  factor  of 
three  larger  that  that  of  S  waves  suggesting  that  the  conventional  assumption  of  dominant  losses 
from  shear  mechanisms  is  not  the  case.  It  should  be  noted  that  these  confining  pressures  are 
small  compared  with  those  for  underground  sources. 

Multicycle  Laboratory  Experiments 

Resonant  Bars.  Tittmann  [1983]  has  taken  laboratory  data  on  the  absorption  of  the  energy  in 
multiple  cycle  oscillations  of  halite  rods.  Both  dome  salt  and  pressed  salt  samples  were  used. 
These  resonant  bar  experiments  measure  the  width  of  the  resonance  peak  for  cyclic  motion  at 
frequencies  from  about  90  to  500  Hz  induced  in  mechanically  loaded  salt  samples.  They  were 
carried  out  for  both  torsion  and  flexure  modes  with  peak  strains  from  10"^  to  10'^.  Pressure 
variation  studies  were  carried  out  using  jacketed  samples  in  pressurized  chambers  allowing 
pressures  up  to  124  MPa.  The  effects  of  humidity  variation  were  also  determined. 
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The  results  for  pressed  salt  indicate  that  attenuation  is  dependent  on  all  parameters  which  were 
varied.  Attenuation  decreases  with  increasing  frequency  and  increases  with  increasing  humidity. 
At  ambient  pressure,  for  strains  below  2x1 0'^  the  attenuation  is  only  weakly  nonlinear  while 
above  this  it  is  strongly  amplitude  dependent.  In  the  low  strain  region  Q's  for  extension  and 
torsion  are  about  500  and  1000,  respectively. 

For  dome  salt  the  results  are  somewhat  different.  Nonlinear  behavior  persists  to  very  low 
strains  (10"^).  For  any  significant  confining  pressures  the  threshold  strain  for  nonlinear 
behavior  increases  to  about  2x10'^  but  the  attenuation  is  then  essentially  independent  of  pressure 
up  to  at  least  68  MPa.  The  attenuation  increases  slightly  with  increasing  frequency;  1/Q  increases 
by  a  factor  of  about  two  from  80  to  480  Hz.  The  high  pressure  behavior  appears  to  stabilize  at 
the  quoted  levels  only  after  a  period  of  adjustment  of  hours  suggesting  that  crack  healing 
strengthens  the  samples  under  pressure.  An  illustration  of  the  strain  dependence  is  shown  in 
Figure  10. 

Tittmann  points  out  that  the  attenuation  in  these  multicycle  experiments  after  the  hundreds  of 
cycles  required,  especially  for  high  strains,  may  reflect  changes  or  damage  in  the  material 
resulting  from  previous  cycles.  Measurements  made  then  may  then  not  correspond  to  behavior 
for  a  single  pulse.  Tittmann  gives  no  data  on  this  subject  but  it  has  been  investigated 
experimentally  by  Bonner  et  al.  [1989]  although  not  with  salt.  The  difference  in  attenuation 
between  uncycled  granite  and  samples  with  10^  cycles,  attributed  to  fatigue  damage,  is  shown  in 
Figure  11.  This  suggests  the  possibility  that  Tittmann's  experimental  strain  dependence  may 
actually  be  the  result  of  damage  from  cycling.  In  any  case,  Tittmann's  attenuation  for  salt  is 
significantly  less  than  that  suggested  by  the  other  experiments  described  in  this  article. 

Summary  Of  Salt  Attenuation 

While  there  is  a  substantial  body  of  data  for  attenuation  of  signals  in  salt,  the  results  are  rather 
diverse.  It  appears  that  attenuation  is  a  function  of  many  factors  including  strain,  frequency, 
humidity,  number  of  cycles,  source  of  salt  samples  and  character  of  experiment.  It  is  difficult,  if 
not  impossible,  to  combine  the  various  experimental  results  into  a  cohesive  pattern,  let  alone  a 
constitutive  relation.  However,  there  is  some  degree  of  consistency  which  we  shall  now  attempt 
to  define. 

It  is  most  striking  that  the  data  from  explosive  sources  roughly  satisfy  simple  cube  root  scaling 
for  peak  velocities,  if  not  for  all  details  of  the  waveforms.  This  means  that  if  the  lengths  and 
times  are  all  scaled  by  the  cube  root  of  the  yield  then,  to  a  remarkable  degree,  all  tamped 
explosive  source  experiments  give  approximately  the  same  pulses  at  any  scaled  range. 

It  is  perhaps  not  so  remarkable  that  this  scaling  should  hold  near  the  explosions  since  the 
initial  pulse  character  is  not  determined  by  attenuation.  Rather  the  scale  of  the  explosively 
generated  cavity  is  fixed  by  the  cube  root  of  ratio  of  the  density  to  the  yield.  Given  the 
propagation  velocity,  this  spatial  scale  will  determine  a  temporal  scale,  both  varying  as  the  cube 
root  of  the  yield. 
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The  more  interesting  result  is  that  the  pulses  continue  to  scale  as  they  propagate  out  into  the 
medium  at  strain  levels,  less  than  10"^,  for  which  there  are  no  gross  changes  in  the  medium.  If 
the  attenuation  suffered  is  intrinsic  to  the  medium,  or  linear,  the  associated  Q  must  be 
independent  of  frequency  in  order  for  the  results  to  obey  cube  root  scaling.  This  is  clear  from  the 
form  shown  in  (1)  for  the  attenuation  operator.  If  scaling  holds,  the  expression  0)r/cQ  must 
scale.  The  combination  cor  has  units  of  velocity  which  scales  while  c  is  constant.  Therefore  Q 
must  scale  but  it  can  only  a  function  of  frequency  since  the  medium  is  assumed  uniform  and 
linear  and  the  only  function  of  frequency  which  scales  is  a  constant.  In  general,  for  either  linear 
or  nonlinear  effects,  the  fact  that  the  experimental  results  cube  root  scale  indicates  that  the 
medium  must  have  no  inherent  scales  of  length  or  time  (in  the  range  of  scales  of  the  experiments) 
so  any  constitutive  relation  must  be  rate  independent. 

Salmon,  Cowboy  and  small  scale  laboratory  velocity  pulses  all  scale  well  but  they  also  all 
indicate  that  the  effective  Q  extracted  is  strongly  frequency  dependent.  To  illustrate  this  clearly, 
consider  Figure  12  which  shows  unsealed  estimates  of  Q  as  a  function  of  frequency  for  pairs  of 
records  from  Salmon  and  Cowboy  at  approximately  the  same  peak  strain.  Note  that  the  two 
functions  of  frequency  are  quite  distinct  but  they  are  quite  similar  when  scaled  relative  to  their 
respective  comer  frequencies.  The  comer  frequency,  of  course,  scales  with  the  cube  root  of  the 
yield  which  says  that  the  scaled  effective  Q's  are  nearly  the  same.  'Hiis  has  been  observed  both 
by  Tmlio  [private  communication]  and  by  Wortman  and  McCartor  [1989].  This  clearly  says  that 
the  attenuation  is  not  just  a  function  of  the  medium  but  it  must  depend  on  amplitude  or  shape  of 
the  pulse.  In  other  words,  the  attenuation  must  be  nonlinear,  at  strains  above  10'^.  Since  the 
behavior  is  nonlinear,  there  is  no  benefit  in  using  a  Q  description.  In  fact,  the  use  of  Q  often 
serves  to  confuse  the  fundamental  problem  of  finding  a  physically  meaningful  constitutive 
relation  at  moderate  strains. 

In  spite  of  the  fact  that  moderate  strain  attenuation  is  almost  certainly  nonlinear,  it  is  possible 
to  use  effective  Q  estimates  to  combine  data  by  taking  the  effective  Q  at  the  dominant  frequency. 
This  effective  Q  gives  a  measure  of  the  magnitude  of  attenuation.  It  is  much  more  difficult  to 
combine  attenuation  information  from  explosive  pulses  and  multicycle  experiments  since  the 
effective  Q  will  generally  be  a  function  of  the  details  of  the  experiment.  If  a  proper  constitutive 
relation  were  known,  a  comparison  could  be  made.  However,  no  such  relation  is  known. 
Ignoring  this  substantial  problem,  the  data  from  all  the  experiments  discussed  in  the  text  of  this 
paper  can  be  expressed  as  effective  Q  and  compared  as  a  function  of  peak  strain.  Figure  13  gives 
all  these  data  on  a  single  plot.  The  strain  range  available  is  from  about  10'^  to  10"^. 

With  a  single  exception,  the  trend  of  this  collection  of  data  is  for  a  decrease  in  attenuation  with 
decreasing  peak  strains.  The  Salmon,  higher  strain  Cowboy  and  small  scale  explosion  data  give 
Q's  of  the  order  ten.  The  Cowboy  and  Sterling  data  suggest  that  Q's  are  well  over  one  hundred 
by  strains  of  10'^.  The  resonant  bar  results  (which  seem  to  consistently  show  lower  attenuation 
than  other  experiments)  have  Q's  of  several  hundred  but  show  an  increase  as  strains  approach 
10'^.  The  New  England  Research  ultrasonic  pulse  experiments  are  the  exception  as  they  give  a 
Q  independent  of  strain.  However,  these  experiments  were  carried  out  with  no  confinement. 
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Furthermore,  recent  work  on  other  media  [K.  Coyner,  private  communication]  suggests  that  this 
experimental  technique  may  be  strongly  influenced  by  scattering  from  normally  existing 
inhomogeneities  in  the  samples.  This  scattering  from  structure  comparable  with  the  ultrasonic 
wavelength  may  depend  strongly  on  frequency  so  the  attenuation  results  would  not  be  relevant  to 
the  lower  frequency  pulse  propagation  problem. 

Taken  as  a  whole,  the  salt  data  indicate  that  attenuation  is  strong  and  nonlinear  for  strains 
greater  that  10"^.  For  smaller  strains  the  attenuation  appears  to  become  small  and  linear.  The 
peak  strain  of  10"^  occurs  at  a  scaled  range  of  approximately  3  km/kt^/^.  Murphy  [1978] 
indicates  that  an  elastic  radius  for  salt  is  less  than  500  m/kt^/^  based  on  analysis  of  stabilization 
of  RDF's  of  nuclear  explosions.  This  contrasts  with  the  fact  that  other  measurements  quoted  in 
the  current  article  indicate  that  the  effective  Q  is  less  than  10  and  apparently  nonlinear.  When  a 
source  function  is  defined  at  this  small  radius  the  subsequent  propagation  out  to  sufficiently  large 
ranges,  where  the  totality  of  the  data  indicates  the  behavior  is  linear,  will  have  an  effect  which  is 
not  currently  understood. 

There  are  no  completely  satisfactory  descriptions  of  the  nonlinear  behavior  of  salt  in  the 
moderate  strain  regime.  Without  a  physical  constitutive  relation,  it  is  difficult  to  remove  the 
uncertainty  in  the  effective  seismic  source  function.  There  are  two  conjectures  which  have  been 
put  forth  to  account  for  Salmon  data. 

Rimer  and  Cherry  [1983]  have  noted  that  a  reasonable  fit  to  the  Salmon  data,  including 
attenuation  and  precursor,  can  be  obtained  by  use  of  a  constitutive  relation  combining  a  limiting 
yield  strength  with  quadratic  work  hardening  and  softening.  Cherry  and  Rimer  [1980]  find  that 
the  same  parameters  also  provide  a  reasonable  description  of  other  salt  data.  The  work  hardening 
model  provides  a  variable  yield  strength  Y  which  is  given  by 

Y  =  Yo(l+eiE-e2E2)<YLta  (3) 

There  are  only  three  free  parameters  since  the  limiting  yield  strength,  YLini’  is  constrained  by 
other  experiments.  In  this  model  the  initial  yield  strength,  Yq,  is  rather  low.  The  harding  and 
softening  parameters  ej  and  e2  as  well  as  Yq  are  used  to  fit  the  data.  The  yield  strength  initially 
increases  as  inelastic  energy,  E,  is  absorbed  by  the  shear  failure  in  a  manner  quadratic  in  this 
inelastic  energy.  The  yield  strength  is  a  measure  of  the  maximum  potential  energy  which  can  be 
held  in  shear  in  the  medium.  When  addition  shear  stress  is  applied  to  a  medium  at  its  yield  limit, 
the  work  goes  into  inelastic  energy  which  then  is  taken  to  alter  the  yield  limit  in  this  model.  As 
indicted  by  Figure  14,  there  are  some  experimental  data  to  suggest  that  work  hardening  and 
yielding  in  salt  does  occur  [Boresi  and  Deere,1963].  However,  Glenn  [1989]  has  noted  that 
while  the  Salmon  data  can  be  largely  accounted  for  by  a  constitutive  model  for  which  the  salt  first 
hardens  and  then  softens  greatly,  the  required  parameters  are  inconsistent  with  independently 
measured  material  properties. 
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McCartor  and  Wortman  [1988,1990]  have  proposed  another  nonlinear  model  for  partial  shear 
failure  which  is  designed  to  account  for  the  Salmon  attenuation  and  precursor.  In  this  model  the 
Lamd  shear  modulus,  p.,  is  permanently  and  instantly  reduced  by  80%  in  any  material  element 
once  the  compressional  strain  level  exceeds  10‘^.  The  other  Lamd  modulus,  X,  is  held  fixed. 
This  has  the  effect  of  reducing  the  compressional  speed  (which  is  proportional  to  +2p)l/2)  of 
the  main  part  of  the  pulse  by  about  20%  relative  to  the  elastic  speed  seen  in  the  precursor.  In 
addition  to  producing  a  precursor  which  propagates  at  the  observed  speed,  this  also  produces 
attenuation  due  to  the  loss  of  shear  energy  which  is  proportional  to  m,  giving  an  effective  Q  of 
about  13  for  peak  strains  well  in  excess  of  10'^.  Note  that  for  small  strains  below  this  threshold 
there  will  be  no  loss.  This  value  of  Q  is  much  less  than  that  expected  for  very  small  strains  but  it 
is  still  more  than  the  5  to  10  seen  for  Salmon  attenuation.  It  will  produce  a  sharply  changing 
effective  Q  at  a  threshold  strain  in  the  manner  suggested  by  the  Cowboy  data.  Furthermore, 
since  the  partial  shear  failure  threshold  is  a  function  of  the  strain  level  only,  the  resulting 
constitutive  model  will  preserve  cube  root  scaling.  McCartor  and  Wortman  [1990]  find  that  this 
mechanism  is  not  adequate  to  account  for  all  aspects  of  the  waveform  unless  some  addition  linear 
attenuation  is  added.  Still  their  calculations,  and  the  model  of  Rimer  and  Cherry  [1983]  (which 
also  is  consistent  with  cube  root  scaling),  strongly  suggest  that  shear  failure  plays  a  strong  role  in 
the  nonlinear  behavior  of  salt  at  moderate  strains.  While  these  models  hint  at  the  character  which 
is  required  for  a  robust  constitutive  relation  for  salt  at  moderate  strains,  the  issue  is  clearly  not 
resolved. 
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Fig.  1.  Salmon  particle  velocity  records  at  166, 225, 276, 318, 402  and  660  m,  in  that 
order  from  left  to  right  [McCartor  and  Wortman,  1989]. 
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Fig.  2.  Ferret's  [1967]  Salmon  peak  radial  particle  velocities  versus  range  from  18 
sensors. 
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Fig.  4.  Salmon  Q's  estimated  from  adjacent  sensor  pairs  between  166  and  660  m 
[McCartor  and  Wortman,  1989]. 
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Fig.  5.  Salmon  Q  estimates  with  range  [Gupta  and  McLaughlin,  1989].  Vertical  bars 
show  a  standard  deviation.  Horizontal  bars  indicate  the  separation  of  the  sensors  used 
for  that  estimate. 
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Fig.  6.  Cowboy  peak  particle  velocity  data  versus  scaled  range  [Murphey.i 


Fig.  7.  Cowboy  Q  estimates  with  range  from  Wortman  and  McCanor  [1989]  and 
MinstCT  and  Day  [1986].  Vertical  bars  show  the  variation  of  Q  over  an  order  of 
magnitude  in  frequency  about  the  dominant  frequency.  Horizontal  bars  indicate  the 
separation  of  the  sensors  used, 
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Fig.  8.  Salt  peak  panicle  velocity  data  versus  scaled  range  [Larson,  1982] 
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Fig.  10.  Resonant  bar  Q  estimates  versus  strain  for  three  pressures  for  both  flexure 
(solid)  and  torsion  (open)  [Tittmann,  1983]. 
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Fig.  1 1.  Variation  of  attenuation  with  strain  for  cycling  in  Sierra  white  granite  [Bonner 
et  al.,1989]. 
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Fig.  12.  Unsealed  Salmon  and  Cowboy  Q  estimates  as  a  function  of  frequency 
[Wortman  and  McCartor,  1989]. 
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1983]).  See  Fig.  10  for  resonant  bar  parameters. 
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